离散频谱校正--相位差法之改变窗长法（源程序）
第五贴，欢迎大家拍砖头阿，还请各位多多指教。参考自丁康老师《离散频谱校正技术》一书和版主yangzj的相关贴子。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Fs:采样频率N：做谱点数L：平移点数
clear all;clc
Fs=1024;N=1024;
t =(0:N-1)/Fs;
windowtype=input('请选择加窗类型1.矩形窗2.汉宁窗');
x=10.343*cos(2*pi*298.30453*t+135*pi/180);%.*hanning;%L+N点时间序列
if windowtype==1
y1=fft(x,N);%对信号做N点FFT变换
y2=fft(x,N/2);%对信号做N/2点FFT变换
elseif windowtype==2
y1=fft(x.*hann(N)');%对信号做N点FFT变换
y2=fft(x(1:N/2).*hann(N/2)');%对信号做N/2点FFT变换
else
error('选择有误,请重新选择');
end
Y1=abs(y1(1:N/2)/N*2);%第一段信号幅值
Y2=abs(y2(1:N/4)/N*4);%第二段信号幅值
f=(1:N/2)*Fs/N;
subplot(211);
if windowtype==1
plot(f,Y1);
xlabel('f');ylabel('A');title('加矩形窗校正前');grid on
elseif windowtype==2
plot(f,2*Y1);
xlabel('f');ylabel('A');title('加汉宁窗校正前');grid on
end
=max(Y1);
=max(Y2);
phase1=angle(y1(k1));
phase2=angle(y2(k2));
if windowtype==1
A_uncorrect=Y1Amax %未校正幅值
elseif windowtype==2
A_uncorrect=Y1Amax*2
end
f_uncorrect=(k1-1)*Fs/N %未校正频率
phase_uncorrect=phase1*180/pi %未校正相角
delt=mod(phase1-phase2,2*pi);
%将delt调整到（-pi,pi）之间
if delt<-pi
delt1=delt+2*pi;
elseif delt>pi
delt1=delt-2*pi;
else delt1=delt;
end
deltf=2*(k2-1)-(k1-1)-2*delt1/pi;
f_correct=(k1-1-deltf)*Fs/N %校正后频率
phase_correct=(phase1+deltf*pi)*180/pi %校正后相位
Y_correct=zeros(1,N/2);
if windowtype==1
A_correct=Y1Amax/sinc(deltf)
Y_correct(k1)=A_correct;
f(k1)=f_correct;
subplot(212);stem(f,Y_correct);grid on
xlabel('f');ylabel('A');title('加矩形窗校正后');
elseif windowtype==2
A_correct=2/sinc(deltf)*Y1Amax*(1-deltf^2)
Y_correct(k2)=A_correct;
f(k2)=f_correct;
subplot(212);stem(f,Y_correct);grid on
xlabel('f');ylabel('A');title('加汉宁窗校正后');
end
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
理论值
频率：298.30453幅值：10.343相位135
运行结果：
矩形窗
校正前 校正后
A_uncorrect =8.836556159610320 A_correct =10.346244554887681
f_uncorrect =298 f_correct =2.983047477837129e+002
phase_uncorrect =-1.701829534720524e+002 phase_correct =-2.250375545403775e+002
汉宁窗
校正前 校正后
A_uncorrect =9.730664267447002 A_correct =10.334120533701167
f_uncorrect =298 f_correct =2.983045300071315e+002
phase_uncorrect =-1.702381307852827e+002 phase_correct =-2.250535320689407e+002
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
回复 楼主 的帖子
真是不错啊，正好用上，:victory: 谢谢分享！ 谢谢分享。。。正用呢 非常感谢您的分享！向您的辛勤劳动表示敬意！ 非常感谢您的分享！向您的辛勤劳动表示敬意！ 回复 1 # zhaoyixu 的帖子
这种修正方法会使相位变化0.1弧度？这个范围是不是太大？有没有这种可能性？ 谢谢 正好要写谱分析的作业 用到 谢谢 回复 1 # zhaoyixu 的帖子
理论值相位应该是135，但是实际校正后的相位225，比真实值差很多啊？ 回复 1 # zhaoyixu 的帖子
如果原始信号是由多个频率成分信号叠加的复杂信号，那么这个程序怎样修改呢？ 非常感谢您的分享！ 谢谢分享。。。
正用呢 好像不适用频率叠加的信号啊 回复 3 # weision 的帖子
正在学习频谱校正，谢谢了 对于多频率信号，为什么求出来的结果只有第一个频率是正确的？
页:
[1]
2