求教：apFFT频谱校正问题！

whoru

我想写一个通用的apFFT时移相位差法校正程序，也就是对任意的信号，可以采用apFFT作频谱校正，主要参考论坛里面贴出的程序，可是得到的结果非常不理想，请问高手是什么问题？

function [fff,AA,PH,N]=apFFT_Correct(s,Fs,CNum)
% 采用时移相位差法校正全相位FFT谱,需要3*N-1个样本点
% close all;
% clear all;
% clc;
N=fix((length(s)+1)/3);
s=s(1:3*N-1);

win=hann(N)';
win1=hann(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);
s1=s1-ones(size(s1)).*mean(s1);
y1=s1.*win2;
y1a=y1(N:end)+[0,y1(1:N-1)];
Out1=fft(y1a,N);
s2=s(1+N:3*N-1);
s2=s2-ones(size(s2)).*mean(s2);
y2=s2.*win2;
y2a=y2(N:end)+[0,y2(1:N-1)];
Out2=fft(y2a,N);
a1=abs(Out1);
p1=rem(atan(imag(Out1)./real(Out1)),2*pi);
a2=abs(Out2);
p2=rem(atan(imag(Out2)./real(Out2)),2*pi);
g=rem((p2-p1)/pi/2,1);
h=2*pi*g.*(1-g.*g)./sin(pi*g);

aa1=abs((h.^2).*a2)/2;
L=fix(N/2)+1;
f=(0:L)*Fs/N;
figure;
subplot(1,2,1)
plot(f,a1(1:L+1));
subplot(1,2,2)
plot(f,p1(1:L+1));
figure;
subplot(1,2,1)
plot(f,a2(1:L+1));
subplot(1,2,2)
plot(f,p2(1:L+1));
for i=1:CNum
    [A,Ind(i)]=max(a1(1:L));
    f0=f(Ind(i));
    fff(i)=f0+g(Ind(i))*Fs/N;
    AA(i)=aa1(Ind(i));
    dp=rem((Ind(i)-1)*2*pi*fff(i)/Fs,2*pi);
    PH(i)=p1(Ind(i))*180/pi;
    StartP=Ind(i)-2;
    StopP=Ind(i)+2;
    if Ind(i)-2<1
        StartP=1;
    end
    if Ind(i)+2>L
        StopP=L;
    end
    Indx=StartP:StopP;
    a1(Indx)=zeros(1,length(Indx));
end




测试程序：
clc
Fs=4096;
N=8192;
t=-N+1:N*2-1;

f=[4.2,126.4,131.7,258.5,382.8];
A=ones(1,5);
Ph=[40,20,30,70,60];
CNum=length(f);
s=zeros(1,length(t));
for i=1:CNum
    myt=A(i)*cos(2*pi*t*f(i)/Fs+Ph(i)*pi/180);
    s=s+myt;
end

[ff,A,pho,Nn,tt]=apFFT_Correct(s,Fs,CNum)
% 利用传统比值法对频谱校正的结果，注意：这里输出的相位值为起始点的相位值。
% N=0点的相位值应等于(输出相位值+2*pi*N*f1/Fs/pi*180)
for i=1:CNum
    Phh(i)=rem(2*pi*Nn*f(i)/Fs+Ph(i)*pi/180,2*pi)*180/pi;
end
[f;pho;ff;Phh]

理想的频率：
    4.2000  126.4000  131.7000  258.5000  382.8000
理想的相位  
70.0000   20.0000   30.0000   40.0000   60.0000
校正的频率：
  3.9500 126.4000 131.4500  258.5000   382.8000
校正的初相位
  184.0000  308.0000  174.0000   70.0000  276.0000

我怀疑是程序中的变量“g”的计算有问题，但是却不知如何修改，请大牛们指点！！

zhwang554

function [fffa,AA,PH]=apFFT_Corret(k0,F2,F1,N);
clc;clear;clf;close all;
Fs=4096;
N=8192;
t=-N+1:N*2-1;
f=[4.2,126.4,131.7,258.5,382.8];
A=ones(1,5);
Ph=[40,20,30,70,60];
CNum=length(f);
s=zeros(1,length(t));
for i=1:CNum
    myt=A(i)*cos(2*pi*t*f(i)/Fs+Ph(i)*pi/180);
    s=s+myt;
end
N=fix((length(s)+1)/3);
s=s(1:3*N-1);
win=hann(N)';
win1=hann(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);
s1=s1-ones(size(s1)).*mean(s1);
y1=s1.*win2;
y1a=y1(N:end)+[0,y1(1:N-1)];
F1=fft(y1a,N);
s2=s(1+N:3*N-1);
s2=s2-ones(size(s2)).*mean(s2);
y2=s2.*win2;
y2a=y2(N:end)+[0,y2(1:N-1)];
F2=fft(y2a,N);
a1=abs(F1);
L=fix(N/2)+1;
f=(0:L)*Fs/N;
for i=1:CNum
    [A,Ind(i)]=max(a1(1:L));
    k=Ind(i);
    peak_F2(i)=F2(k);
    peak_F1(i)=F1(k);
    StartP=Ind(i)-2;
    StopP=Ind(i)+2;
    if Ind(i)-2<1
        StartP=1;
    end
    if Ind(i)+2>L
        StopP=L;
    end
    Indx=StartP:StopP;
    a1(Indx)=zeros(1,length(Indx));
end
num=length(Ind);  %%%%%设置所取的“谐波簇”个数
for m=1:num;        
     [fff(m) ,AA(m) ,PH(m) ]= apFFT_Corret(Ind(m),peak_F2(m),peak_F1(m),N);
end
disp('校正的频率')
fm=[fff]*Fs/N
disp('校正的相位')
PH=mod(PH*180/pi,360)
disp('校正的振幅')
AA
function [fff ,AA ,PH ]= apFFT_Corret(k0,F2,F1,N)   
dph=angle(F2)-angle(F1);
dph=mod(dph,2*pi);
     if dph>pi
        dph=dph-2*pi;
      elseif dph<-pi
        dph=dph+2*pi;
     end   
     dk=dph/pi/2;
     g=dk;
      h=2*pi*g.*(1-g.*g)./sin(pi*g);
      AA=abs((h.^2).*F1)/2;
      fff=(k0+dk-1);
      PH=angle(F1);

运行结果:
校正的频率fm               258.5000  126.4000  131.7000    4.2000  382.8000
校正的相位PH              70.0000    20.0000   30.0000   40.0000    60.0000
校正的振幅AA                1.0000     1.0000     1.0001      1.0001      1.0001

[ 本帖最后由 zhwang554 于 2009-6-29 18:54 编辑 ]

dsfire

王老师，您好，您上面的程序中，若f的个数是未知的，那CNum该怎么求呢？

dsfire

您好，这个程序里面进行了来回调用，难道我理解错了？

zhwang554

程序的搜索语句部分是楼主写的, 我只是在上面加了判断语句

可设定一门限,搜索振幅大於该值的峰值

[ 本帖最后由 zhwang554 于 2010-1-22 17:15 编辑 ]

dsfire

明白了！谢谢！再请问一下，执行下面的语句，程序不就来回调用了么？
for m=1:num;        
     [fff(m) ,AA(m) ,PH(m) ]= apFFT_Corret(Ind(m),peak_F2(m),peak_F1(m),N);
end

zhwang554

对呀,把找到的峰值,一个一个调用apFFT作校正

dsfire

那来回调用，程序不就死循环了么？上面的程序我直接运行，程序进入死循环了。
我学习matlab不多，对function认识不够，还望老师指教！谢谢老师了！

zhwang554

num=length(Ind);  %设置所取的“谐波簇”个数
m=1:num;     
到num就完了
我刚调用上程序, 运行没问题呀, 你一定改动了,若没有把握就不要改

[ 本帖最后由 zhwang554 于 2010-1-22 20:16 编辑 ]

dsfire

谢谢老师的回复！
我是直接复制粘贴到一个文件里运行的，是不是这出错了！

dsfire

老师，能把复制好的文件发给我一份么？我这里运行就死循环了。zisehuoxia@126.com，谢谢了！:@)

zhwang554

就是沙发中的程序, 运行没问题

dsfire

谢谢老师，搞定了！复制错了！:@)

[ 本帖最后由 dsfire 于 2010-1-23 14:41 编辑 ]

huan0918

但是我把输入频率该大点就出问题了
f=[5001,4006.4,1231.7,2500.5,1582.8];幅度【1，2，3，4，5】
结果是
校正的频率

fm =
  1.0e+003 *
    1.5828    1.5955    1.2317    0.0896    0.9050
校正的相位
PH =
   60.0000  290.0000   30.0000  340.0000   40.0000
校正的振幅
AA =
    5.0000    4.0000    3.0000    2.0000    1.0000
请教老师是哪里出了问题

huan0918

你好，我想问下，要是要校正的频率不知道怎么办，而且是多频的

比方说要分析的信号可能是0-10KHZ的，采样频率是40KHZ

是不是这个方法就不行了
给个建议