学写程序之三(源代码)－－FFT+FT连续细化傅立叶变换分析校正法
第三贴，欢迎大家拍砖头阿，还请各位多多指教。参考自丁康老师《离散频谱校正技术》一书和版主yangzj的相关贴子。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear all;clc
Fs=1024;N=1024;
t =(0:N-1)/Fs;
tt=0:N-1;
hanning=0.5-0.5*cos(2*pi*tt/N);
windowtype=input('请选择加窗类型1.矩形窗2.汉宁窗');
if windowtype==1
x=4.2366*cos(2*pi*63.2*t+23.8*pi/180);%时域仿真函数
elseif windowtype==2
x=4.2366*cos(2*pi*63.2*t+23.8*pi/180).*hanning;
else
error('选择有误,请重新选择');
end
y = fft(x);
Y=y(1:N/2+1)/N*2;
f=(0:N/2)*Fs/N;
subplot(211);plot(f,abs(Y(1:N/2+1)));grid on
A=abs(Y);
=max(A);
if windowtype==1
Amax_uncorrect=Amax%校正前幅值
elseif windowtype==2
Amax_uncorrect=Amax*2
end
phmax_uncorrect=angle(Y(k))*180/pi%校正前相位
f_uncorrect=k-1%校正前频率
L=80;%要细化的点数
deltf=((k+1-1)*Fs/N-(k-1-1)*Fs/N)/L;%确定细化后频率分辨率
YY=zeros(1,N);
ff=(k-1-1)*Fs/N;
for i=1:L
for ii=1:N
YY(i)=YY(i)+x(ii)*exp(-j*2*pi*ff*(ii-1)/Fs);
end
ff=ff+deltf;
end
if windowtype==1;
A_correct=max(abs(YY)/N*2)%校正后幅值
=max(YY);
f_correct=(k-1-1)*Fs/N+(YYk-1)*deltf%校正后频率
phmax_correct=angle(YYmax)*180/pi%校正后相位
Y(k)=A_correct;f(k)=f_correct;
subplot(212);plot(f,abs(Y(1:N/2+1)));grid on
else windowtype==2;
A_correct=2*max(abs(YY)/N*2)%校正后幅值
=max(YY);
f_correct=(k-1-1)*Fs/N+(YYk-1)*deltf%校正后频率
phmax_correct=angle(YYmax)*180/pi%校正后相位
Y(k)=A_correct/2;f(k)=f_correct;%幅值除以2,保证向里的数值为原来数值
subplot(212);plot(f,2*abs(Y(1:N/2+1)));grid on%2倍为幅值恢复系数
end
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
运行结果：
理论值
幅值：4.2366 频率：63.2 相位：23.8度
加矩形窗
校正前
Amax_uncorrect =3.962471854470429
phmax_uncorrect =59.672485657326710
f_uncorrect =63
校正后
A_correct =4.235275753915142
f_correct =63.200000000000003
phmax_correct =23.660397117658885
加汉宁窗
校正前
Amax_uncorrect = 4.128431340024557
phmax_uncorrect =59.800004751788116
f_uncorrect =63
校正后
A_correct =4.236600313219102
f_correct =63.200000000000003
phmax_correct =23.800007456676489 校正后，还是挺准的啊。为何没有人回复呢？现在大家一般用哪种方法进行校正呢？ 学习了多谢
页:
[1]