|
回复 楼上 czk108 的帖子
p22(k)=p2(ff+1) 是 弧度 p22(k)=p2(ff+1)*180/pi 是 度. 程序中ph1=60 度, 所以用 度.
频率取整数时误差为0.5还正不好介决, 到不是判断语句的范围取的不对, 无噪声p1-p2应是零,有噪时p1-p2可能是正一点点,也可能是负一点点, 当负一点点时输出就多了1
对fft/apfft法, 可不用判断语句, 因为fft的相位谱不是水平相位特性,频偏>0.5或<0.5相位不同, 自行判断了. 如下程序
close all;clc;clear all;
T=1;k=1;
while T>0;
T=T-1;
N=1024;
f1=155.72;
A1=1;
ph1=60;
w=2*pi;
t=-N+1:N-1;
y=A1*cos(1*(2*pi*t*f1/N+ph1*pi/180));
y=awgn(y,15);
y1 = y(N:2*N-1);%后N个输入数据
win = hanning(N)';;
win1 = win/sum(win);%窗归1
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);%FFT振幅谱
p1 = mod(phase(y11_fft)*180/pi,360);;%FFT相位谱
y2 = y(1:2*N-1);%2N-1个输入数据
win = hanning(N)';;
winn = conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);%窗归1
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的apFFT输入数据
y2_fft = fft(y222,N);;
a2 = abs(y2_fft);%apFFT振幅谱
p2=mod( phase(y2_fft)*180/pi,360);%apFFT相位谱
ee=(p1-p2)/180/(1-1/N);%频率偏离校正值
aa=(a1.^2)./a2*2;%振幅校正值
rr=round(f1);
p22(k)=p2(rr+1);
aaa(k)=aa(rr+1);
f22(k)=rr+ee(rr+1);
k=k+1;
end
disp('phase correction')
p22=sqrt(mean(p22.^2))
disp('frequency correction')
f22=sqrt(mean(f22.^2))
disp('amplitude correction')
aaa=sqrt(mean(aaa.^2))
你试用上程序,有否问题
ff=round(f1) 中round语句是4舍5入,ff+1点是振幅峰值点.+1是因为Matlab数值从1始,
当频率分辨率不为1时, 如采样频率为fs时, ff=round(f1*N/fs);
[ 本帖最后由 zhwang554 于 2010-1-15 03:22 编辑 ] |
|