离散频谱校正--相位差法之改变窗长法（源程序）

zhaoyixu

第五贴，欢迎大家拍砖头阿，还请各位多多指教。
参考自丁康老师《离散频谱校正技术》一书和版主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
[Y1Amax,k1]=max(Y1);
[Y2Amax,k2]=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
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

guiguihong

真是不错啊，正好用上，:victory:

weision

谢谢分享！

robustcococole

谢谢分享。。。
正用呢

wzwry

非常感谢您的分享！向您的辛勤劳动表示敬意！

wzwry

非常感谢您的分享！向您的辛勤劳动表示敬意！

qiuyun0214

回复 1 # zhaoyixu 的帖子

这种修正方法会使相位变化0.1弧度？这个范围是不是太大？有没有这种可能性？

wangzheyou12

谢谢 正好要写谱分析的作业 用到 谢谢

shajue

回复 1 # zhaoyixu 的帖子

理论值相位应该是135，但是实际校正后的相位225，比真实值差很多啊？

shajue

回复 1 # zhaoyixu 的帖子

如果原始信号是由多个频率成分信号叠加的复杂信号，那么这个程序怎样修改呢？

Doer

非常感谢您的分享！

dongyuanzhixing

谢谢分享。。。
正用呢

yaotaiyang

好像不适用频率叠加的信号啊

caiyuan008

回复 3 # weision 的帖子

正在学习频谱校正，谢谢了

ab77977

对于多频率信号，为什么求出来的结果只有第一个频率是正确的？