求助!频谱校正法

flyinghua

烦劳各位大侠帮我修改个简单的程序;谢谢
%%%不补偿，根据离散频谱时移相位差校正法(多个频率),取点从1:2N.第一段取1:N;第二段取N+1:2N
clc;clear;close all;
N=1024;
w=2*pi/N;
m=1:6;
f(m)=[15.2 50 150 156.5 250 350]         %%%实际频率
n=1:2*N;
a1=0.08*310;a2=310;a3=0.2*310;a4=0.05*310;a5=0.14*310;a6=0.12*310;
x=a1*sin(w*f(1)*n+10*pi/180)+a2*sin(w*f(2)*n+5*pi/180)+a3*sin(w*f(3)*n+20*pi/180)+a4*sin(w*f(4)*n+150*pi/180)+a5*sin(w*f(5)*n+120*pi/180)+a6*sin(w*f(6)*n-150*pi/180);  %%%采样序列  2n个点
x1=x(1:N);   %%%序列1
win=hanning(N)';
xin1=x1.*win;  %%%加窗
     Faw1=fft(xin1);  %%%是除窗还是除N
    phase_ap1=angle(Faw1);
    subplot(2,1,1)
    stem(abs(Faw1),'.')
   
x2=x(N+1:2*N);  %%%序列2
xin2=x2.*win;     %%%加窗
    Faw2=fft(xin2);
    phase_ap2=angle(Faw2);
    subplot(2,1,2)
    stem(abs(Faw2),'.')  %%%纵坐标是？？
%%%
kk=1:N
  delta_phase(kk)=phase_ap2(kk)-phase_ap1(kk)-2*pi*kk;
   ph(kk)=mod(delta_phase(kk),2*pi);
  if ph(kk)>pi
      ph(kk)=ph(kk)-2*pi;
  else if ph(kk)<-pi
      ph(kk)=ph(kk)+2*pi;
  end
  delta(kk)=ph(kk)/(2*pi);     %%%归一化的频谱校正量  
fre(kk)=kk+delta(kk);  %%% 估计频率
Ri=real(Faw1);Ii=imag(Faw1)
phase_ap=arctan(Ii/Ri)+delta_phase(kk)  %%% 估计相位
Ai=faw1/W(delta)   %%% 估计幅值  （ W(delta) 是窗函数的频谱模函数）
disp('真实频率')
f(m)
disp('估计频率')  %%%估计频率有1024个？？？
fre(kk)
disp('估计相位')
phase_ap
disp('估计幅值')   %%%
Ai

yangzj

做校正的时候不需要对所有的谱线都做校正，而是在有频率的地方才做校正，这些频率用寻峰得到

songzy41

版主已在http://forum.vibunion.com/forum/thread-41740-1-1.html中给出了重心法的修正程序，楼主可参考一下修改自已程序。

[ 本帖最后由 zhlong 于 2007-6-14 16:25 编辑 ]