APFFT之后怎么会偏差这么大？

鹿鸣

N=1024;
w=2*pi;
t=-N+1:N-1;
y=220*exp(j*(w*t*50.5/N+60*pi/180))+0.6*exp(j*(w*t*101/N+45*pi/180))+25*exp(j*(w*t*151.5/N+80*pi/180))+...
+0.5*exp(j*(w*t*202/N+42*pi/180))+7*exp(j*(w*t*252.5/N+65*pi/180))+0.4*exp(j*(w*t*303/N+35*pi/180))+4*exp(j*(w*t*353.5/N+120*pi/180))+0.3*exp(j*(w*t*404/N+134*pi/180))+2*exp(j*(w*t*454.5/N+62*pi/180));
y1 = y(N:end);
win =  hanning(N)';;
win1 = win/sum(win);
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);
p1 = mod(phase(y11_fft)*180/pi,360);
y2 = y(1:2*N-1);
win =  hanning(N)';;
winn =  conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的apFFT输入数据
y2_fft = fft(y222,N);
a2 = abs(y2_fft);
p2=mod( phase(y2_fft)*180/pi,360);
      ee=(p1-p2)/180/(1-1/N);
      aa=(a1.^2)./a2;
subplot(4,1,1),stem(a2,'.');title('apFFT振幅谱');ylim([0,1]);xlim([0 N/2]);grid
subplot(4,1,2),stem(p2,'.');title('初相位校正谱');ylim([0,400]);xlim([0 N/2]);grid
subplot(4,1,3);stem(ee,'.');title('频率校正谱');ylim([-1,1]);xlim([0 N/2]);grid
subplot(4,1,4);stem(aa,'.');title('振幅校正谱');ylim([0,1.5]);xlim([0 N/2]);grid
      disp('相位校正值')
相位校正值
     p2(50:100:450)
ans =
   60.0000   80.0000   65.0000  120.0000   62.0000
      disp('频率校正值')
频率校正值
    ee(50:100:450)+(49:100:450)
ans =
   50.5000  150.4990  250.4980  348.4951  450.4959
      disp('振幅校正值')
振幅校正值
      aa(50:100:450)
ans =
  220.0001   24.9944    6.9969    3.9977    1.9995

zhwang554

      你没有在频谱的各成份的振幅峰值处取值, 特别是频率校正的整数值必须是峰值频率

      fft/apfft校正法的相位校正值,振幅校正值在峰值附近都有一个平台, 所以在峰值附近取值都可以,但在峰值处取值最好


相位校正值
     p2(50:100:450)
频率校正值
    ee(50:100:450)+(49:100:450)
振幅校正值
   aa(50:100:450)

改为

相位校正值
     p2(51:101:455)
频率校正值
    ee(51:101:455)+(50:101:455)
振幅校正值
      aa(51:101:455)

鹿鸣

非常感谢王老师的指点，我最近也在研读你写的书，非常受用！

鹿鸣

还有一个问题，下面的这个段程序，为什么就是倒数第二个数据相差的很大，其他数据都是准确的，希望指点迷津。
clc;clear;clf;close all;
f0=50.2;
Fs=50*128;
N=1024;%2048/256=8
%n=[0:1:N-1];
n=-N+1:N*2-1;
A=[220,0.6,25,0.5,7,0.4,4,0.3,2];
ph=[60,45,80,42,65,35,120,134,62];
f=[50.5,101,151.5,202,252.5,303,353.5,404,454.5];
s=zeros(1,length(n));
for i=1:9
myt=A(i)*cos(2*pi*f(i)*n/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:9
    [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);

校正的相位
PH=mod(PH*180/pi,360)
PH =
  Columns 1 through 8
   60.0000   80.0000   65.0000  120.0000  180.0000   62.0000   45.0004   42.0006
  Column 9
   35.0005
disp('校正的振幅')
校正的振幅
AA
AA =
  Columns 1 through 8
  220.0035   25.0037    7.0029    4.0020    2.8821    2.0004    0.6000    0.5001
  Column 9
    0.4002
>> PH =
  Columns 1 through 8
   60.0000   80.0000   65.0000  120.0000  180.0000   62.0000   45.0004   42.0006
  Column 9
   35.0005
只有一个数据有问题，而且都是倒数第二个，很奇怪也很有规律。

鹿鸣

校正的频率
fm=[fff]*Fs/N
fm =
  Columns 1 through 8
   50.5000  151.5000  252.5000  353.5000    3.1250  454.5000  101.0000  202.0000
  Column 9
  303.0000
忘记贴频率的了。

zhwang554

去除语句
s1=s1-ones(size(s1)).*mean(s1);
s2=s2-ones(size(s2)).*mean(s2);
即可

鹿鸣

谢谢王老师的指点迷津，让我茅塞顿开。

鹿鸣

王老师你好，能不能解释下为什么去掉
s1=s1-ones(size(s1)).*mean(s1);
s2=s2-ones(size(s2)).*mean(s2);
就能得到正确的结果，我没有想明白，谢谢。

zhwang554

　　原地板贴子中的结果不是倒数笫2个错了, 而是漏检了,
　　你读出地板程序中的Ind值(笫52行后加Ind=Ind),
　　Ind =9 25 41 58 1 74 17 33 9
　　这是信号在N=1024,Fs=50*128时的apfft频谱图, 对照原信号,出错的频率为Ind=1, 即等效直流份量(N=1024,Fs=50*128时的apfft频谱图中的), 也就是说, 信号中的等效直流分量也作为一个峰值被检出了,
　　我用日志<自动搜索峰值的 apfft/apfft 和 fft/apfft 校正程序>中的程序和另一个自动搜索峰值程序分析你的信号,都没有这个问题.对比程序是多了这二个语句.
　　下图a为不加去直流语句的频谱,下图b为加去直流语句的频谱.从图b可见,有一直流份量峰值.
　　这二个语句本耒是为了去掉直流分量不彻底遗留下的, 去掉上面二个语句以后等效直流峰值就没有了,搜索语句就将最小的那个频率成分峰值搜索出耒了.所以对原始信号作fft前不要任意作予处理, 如s1=s1-ones(size(s1)).*mean(s1) 语句在离散和截断情况下, 不一定能彻底去直流.

[img=28,30][/img]
　　还有一个办法是全部地板程序不改动(保留上二语句),将第37句
　　for i=1:9
　　改为 for i=1:10;
　　则多搜索一个峰值, 这样虽直流份量占了一个峰值, 搜索出第10个峰值即原漏检的那个
　　事实上,我们分析一个信号並不知它有几个峰值, 一个实用的自动搜索峰值校正程序应将大於某一阈值的峰值全部搜索出来校正,供使用者分析判断

鹿鸣

这下理解了，多谢王老师仔细的讲解，呵呵。

鹿鸣

我把地板帖子中的汉宁窗改成布莱克曼窗或者布莱克曼哈里斯窗，得到的频率和相位偏差不大，但是幅值却差了好多，是不是可以理解为汉宁窗是最适合用于APFFT的呢。用布莱克曼窗替代汉宁窗之后得到的幅值如下：
校正的振幅
AA
AA =
  Columns 1 through 8
  220.4059   25.4221    7.3396    4.2375    2.0463    0.6045    0.5152    0.4286
  Column 9
    0.3117

zhwang554

　　Blackman等其它窗不能在地板apfft/apfft校正程序中用, 地板程序中的插值公式只针对hanning窗. Blackman窗另有插值公式计算, 效果也不错. 不能把地板帖子中的汉宁窗改成布莱克曼窗或者布莱克曼哈里斯窗耒校正振幅.
　　但fft/apfft校正法中你可将汉宁窗改成布莱克曼窗或者布莱克曼哈里斯窗耒校正振幅.

  
   最好是1.win=hann   win1=hanning, 比2..win=hann  win1=hann 好2个数量级, 比3..win=hanning  win1=hanning 也好2个数量级, 你注意2 win=hann  win1=hann的校正振幅都大於原值, 3 win=hanning  win1=hanning的都小於原值. 所以两者的结合1 win=hann  win1=hanning的精度最好.
1. win=hann   win1=hanning
219.999998338385     24.9999984573797      6.99999894713531     3.9999991582252      
1.99999984416229     0.600008884400799    0.500006695385548    0.400003561427089   
0.300001500976019
2. win=hann  win1=hann
220.003548407911     25.0036535183861      7.00288078906281     4.00200086192488   
2.00039892518265     0.600048742341943    0.500137676994756    0.400243067138395   
0.300101136181152
3. win=hanning  win1=hanning
219.996448321618     24.9963439427626      6.99711829308886     3.99799845593329      
1.99960084281827     0.599969248273551    0.499875765417038   0.399764200986045  
0.299901899235198

鹿鸣

哦，原来这样，我把书再好好看看，受益匪浅啊，非常感谢王老师的指导。