全相位书中附录3遇到出现幅值NaN的问题

无限剑制

论坛中给的校正的判断部分是:
%%%%%%%%%%%%%%%%%%%%%
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;
%%%%%%%%%%%%%%%%%%%%

这里如果g=0会出现NaN的结果
在丁康老师的书里好象叫什么死点的,这个我不太懂，希望各位高手能给点建议！！


下边是程序：
function  [] = apFFT_Correct(N,Fs,f,A,Ph)
    CNum = length(f);
    t=-N+1:N*2-1;          %apfft时移相位差法需要3N-1点
    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
plot(s);
s=s(1:3*N-1);
win=hanning(N)';
win1=hann(N)';
win2=conv(win,win1);   
win2=win2/sum(win2);   %归一化卷积窗
w=pi*2;
s1=s(1:2*N-1);     %取前2N-1点做FFT
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);   %后2N-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;  
      disp('振幅校正前')
      peak_F2(m)
      peak_F1(m)
      [fff(m) ,AA(m) ,PH(m) ]= apFFT_Corret1(Ind(m),peak_F2(m),peak_F1(m));
end
disp('校正的频率')
fm=[fff]*Fs/N
disp('校正的相位')
PH=mod(PH*180/pi,360)
disp('校正的振幅')
AA
function [fff ,AA ,PH ]= apFFT_Corret1(k0,F2,F1)   
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);

下边是测试程序：

%%%%%%%%%%%%%%%%%%校正程序测试
close all;
clear all;
clc;
%%%%%%%%%%%%%%%%%%%%%校正方法选择：0，时移位相位差法；1，fft/apfft法   %%%%%%%%%%%%%
       method = 0;                                                      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs = 200;  %采样频率Hz
N = 256;   %FFT计算的点数,实际采样点数比这个多2N-1/3N-1           
f1  =  [19.3,40,100];    %频率
a1  =  [0.8,1,1];        %幅度
ph1 =  [100,50,120];     %相位
if method
fft_apfft(N,Fs,f1,a1,ph1)
else
apFFT_Correct(N,Fs,f1,a1,ph1);
end

%这里我模拟采样长度不是整树周期，Fs=200, N 固定255点。幅度校正出了个NaN

无限剑制

那个Fs = 200 N=256,上边打错了

zhwang554

f1=100, 在fs=200, N=256时, 折合频率=f1*N/fs=128, 即振幅谱峰值在N=256频谱的中点,不只振幅校正不对,相位也不对.

在其它点, g=0(折合频率是整数)没有关系, 仍能校正振幅


[ 本帖最后由 zhwang554 于 2010-2-26 17:49 编辑 ]

无限剑制

哦，明白了，我这没满足采样定理，信号没学好。
感谢王老师的指点:lol

jinyi7016

能把 apFFT_Correct 的程序分享下么？感谢