whoru 发表于 2009-6-23 16:33

求教：apFFT频谱校正问题！

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

function =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)+;
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)+;
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
    =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=;
A=ones(1,5);
Ph=;
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

=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


理想的频率：
    4.2000126.4000131.7000258.5000382.8000
理想的相位
70.0000   20.0000   30.0000   40.0000   60.0000
校正的频率：
3.9500 126.4000 131.4500258.5000   382.8000
校正的初相位
184.0000308.0000174.0000   70.0000276.0000

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

zhwang554 发表于 2009-6-29 18:51

apfft/apfft校正程序

function =apFFT_Corret(k0,F2,F1,N);
clc;clear;clf;close all;
Fs=4096;
N=8192;
t=-N+1:N*2-1;
f=;
A=ones(1,5);
Ph=;
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)+;
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)+;
F2=fft(y2a,N);
a1=abs(F1);
L=fix(N/2)+1;
f=(0:L)*Fs/N;
for i=1:CNum
    =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;      
   = apFFT_Corret(Ind(m),peak_F2(m),peak_F1(m),N);
end
disp('校正的频率')
fm=*Fs/N
disp('校正的相位')
PH=mod(PH*180/pi,360)
disp('校正的振幅')
AA
function = 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.5000126.4000131.7000    4.2000382.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 发表于 2010-1-21 14:20

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

dsfire 发表于 2010-1-22 15:42

回复 沙发 zhwang554 的帖子

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

zhwang554 发表于 2010-1-22 16:46

回复 板凳 dsfire 的帖子

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

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


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

dsfire 发表于 2010-1-22 17:21

回复 5楼 zhwang554 的帖子

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

zhwang554 发表于 2010-1-22 18:10

回复 楼上 dsfire 的帖子

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

dsfire 发表于 2010-1-22 19:43

回复 7楼 zhwang554 的帖子

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

zhwang554 发表于 2010-1-22 19:51

回复 楼上 dsfire 的帖子

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

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

dsfire 发表于 2010-1-22 20:07

回复 9楼 zhwang554 的帖子

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

dsfire 发表于 2010-1-23 09:24

回复 9楼 zhwang554 的帖子

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

zhwang554 发表于 2010-1-23 12:12

回复 楼上 dsfire 的帖子

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

dsfire 发表于 2010-1-23 13:27

回复 12楼 zhwang554 的帖子

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

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

huan0918 发表于 2010-3-23 22:05

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

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

huan0918 发表于 2010-3-24 18:17

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

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

是不是这个方法就不行了
给个建议
页: [1] 2
查看完整版本: 求教：apFFT频谱校正问题！