将1周期延拓12周期，APFFT分析出错！

dsfire

采样频率是3200Hz，信号频率为50.4Hz。数据1是直接采样12个周期，进行APFFT/APFFT分析，分析结果没有问题；数据2是采样1周期，然后延拓为12周期，分析结果就出问题了，望赐教！
数据1分析结果：
校正的频率
50.400000 100.799993 151.200000 201.599977 252.000000 302.399976 352.800000 403.199984 453.600001        
校正的相位
60.000000 45.004118 79.999996 42.006698 65.000015 119.997622 194.999960 284.989180 61.999874        
校正的幅值
219.999997 0.600092 24.999998 0.500089 7.000001 0.400044 3.999997 0.299981 1.999991
数据2分析结果：
校正的频率
50.000000 100.000000 150.000000 200.000000 250.000000 300.000000 350.000000 400.000000 450.000000        
校正的相位
55.580539 265.888827 67.087927 272.277699 40.707085 213.682722 176.251443 263.722936 10.007477        
校正的幅值
  NaN NaN NaN NaN NaN NaN NaN NaN NaN

运行数据1时将数据2注释掉，运行数据2时，将数据1注释掉。谢谢！
程序如下： (有附件)

clc;clear;clf;
close all;
format short;
Fs=3200;        %采样频率
N=3200/50*4;    %采样长度
t=(-N+1):(N*2-1);
f0=50.4;        %数据频率
ff=[f0*1 f0*2 f0*3 f0*4 f0*5 f0*6 f0*7 f0*8 f0*9];
Ph=[60,  45,  80,  42,  65,  120,  195, 285, 62];
A= [220, 0.6, 25,  0.5, 7,   0.4, 4,   0.3, 2];
CNum=length(A);
s=zeros(1,length(t));
% *******运行数据1的时候将数据2注释掉，运行数据2的时候将数据1注释掉****
%
% 数据1 没有问题，但数据2分析的就出现问题了！频率不是50Hz时(50.4)，数据2分析时仍然是50Hz，相位和幅值也不对，望赐教
% 数据1
% % % % % % % % % % 数据1：直接采样12个周期，其值在s里面% % % % % % % % % % % % % % %
% for i=1:CNum
%     myt=A(i)*cos(2*pi*ff(i)*t/Fs+Ph(i)*pi/180);
%     s=s+myt;
% end
% % % % % % % % % % % % % % %   数据1 构造结束% % % % % % % % % % % % % % % % % % %
% figure(1);plot(s);
% 数据2
% % % % % % % % % % 数据2：先采样一个周期，再延拓为12个周期，其值在s里面 % % % % % % %
%先采样一个周期数据
N=64;
t=0:N-1;
s2=zeros(1,length(t));
for i=1:CNum
    myt=A(i)*cos(2*pi*ff(i)*t/Fs+Ph(i)*pi/180);
    s2=s2+myt;
end
N=3200/50*4;
% 开始延拓为12个周期
for i=1:3*N-1
    if(i <=64)
        s(i) = s2(i);
    else
        s(i) = s(i-64);
    end
end
% figure(2);plot(s);
% % % % % % % % % % %% % % %    数据2 构造结束   % % % % % % % % % % % % % % % % % %

% % % % % % % % % % % % 以下开始对s进行分析

N=fix((length(s)+1)/3);
s=s(1:3*N-1);
win=hann(N)';
win1=hanning(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
s1=s(1:2*N-1);
y1=s1.*win2;
y1a=y1(N:end)+[0,y1(1:N-1)];
F1=fft(y1a,N);

s2=s(1+N:3*N-1);
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;
i=0;
for I=2:length(a1)/2-1;
    k=I;
    if ( (a1(k)>a1(k-1)) & (a1(k)>a1(k+1)) );
        [aaa,bbb,ccc]= apFFT_Corret(k,F2(k),F1(k),N);%%%aPFFT_Correct函数在下面
        if i<9          %只显示前9次谐波
            i=i+1;
            fff(i)=aaa;
            AA(i)=bbb;
            PH(i)=ccc;
            Ind(i)=k;
        end
    end
end
num=length(Ind);  %%%%%设置所取的“谐波簇”个数
% disp('校正的频率')
fm=[fff]*Fs/N;
% disp('校正的相位')
PH=mod(PH*180/pi,360);
disp('校正的频率')
fprintf('%05.6f\t',fm);
fprintf('\n');
disp('校正的相位')
fprintf('%05.6f\t',PH);
fprintf('\n');
disp('校正的幅值')
fprintf('%05.6f\t',AA);
fprintf('\n');

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);

apFFT_apFFT_test.m (2.31 KB, 下载次数: 4)

2010-9-30 10:24 上传

点击文件名下载附件
测试程序

apFFT_Corret.m (277 Bytes, 下载次数: 4)

2010-9-30 10:24 上传

点击文件名下载附件
apFFT_Correct函数

zhwang554

    观察数据1 和数据2 , 数据1在f0=50.4时不是以64为周期的信号, 和数据2不同,所以apfft后数据1结果正确, 数据2不正确

数据1

145.3200  116.0112   85.4938   55.9640   28.7347    3.8960  -19.2539  -41.2885  -61.9095  -79.9657

-94.2809 -104.7330 -112.7784 -120.9503 -131.5509 -145.3389 -161.0627 -176.1630 -188.2142 -196.2100

-200.9300 -204.2408 -207.8642 -212.4107 -217.2039 -220.8300 -221.8919 -219.4248 -212.8167 -201.5215

-185.0025 -163.0901 -136.5003 -107.0217  -77.0545  -48.6724  -22.7923    0.9536   23.5133   45.3489

  65.7890  83.3247    96.7120  106.1295  113.5288  121.8149  133.2849 148.2866   164.9649  180.2815

191.6900  198.4483   201.8150  204.1087  207.2945  211.9703  217.2543 221.4378   222.8261  220.2158

212.8688  200.2771   182.1253  158.5921  130.7331  100.5113   70.2603  41.8129    15.8548   -8.0233

-30.5794  -51.9982  -71.4882   -87.7475 -100.0176 -108.9734 -116.7456 -125.9362 -138.1903 -153.2291

-168.9844 -182.7717 -192.7715 -198.9112 -202.6617 -205.9714 -210.0744 -214.8955 -219.2872 -221.7618

-221.1230 -216.6189 -207.7015 -193.8094 -174.5272 -150.0794 -121.7417  -91.7114  -62.3495  -35.2121

数据2

115.0483    84.5312   55.0607   27.9119    3.1395  -19.9694  -41.9692  -62.5295  -80.4820  -94.6685

-105.0137 -113.0215 -121.2404 -131.9411 -145.8197 -161.5665 -176.6025 -188.5306 -196.4020 -201.0459

-204.3442 -207.9948 -212.5655 -217.3451 -220.9095 -221.8717 -219.2814 -212.5332 -201.0798 -184.3888

-162.3111 -135.5969 -106.0649  -76.1207  -47.8113  -22.0103    1.6832   24.2180   46.0259   66.3990

  83.8172  97.0660   106.3811  113.7627  122.1222  133.7114  148.8050  165.4886  180.7132  191.9741

198.5964  201.8916   204.1888  207.4212  212.1366  217.4138  221.5327  222.8096  220.0584 2 12.5520

199.7871  181.4579   157.7659  129.7977  115.0483  84.5312   55.0607    27.9119    3.1395  -19.9694

-41.9692  -62.5295   -80.4820  -94.6685 -105.0137 -113.0215 -121.2404 -131.9411 -145.8197 -161.5665

-176.6025 -188.5306 -196.4020 -201.0459 -204.3442 -207.9948 -212.5655 -217.3451 -220.9095 -221.8717

-219.2814 -212.5332 -201.0798 -184.3888 -162.3111 -135.5969 -106.0649  -76.1207  47.8113   -22.0103

数据1在f0=50.0时是以64为周期的信号, 和数据2相同,这时apfft后数据1结果正确, 数据2也正确

dsfire

多谢老师解答，十一刚回来，不好意思。对于这种不是整周期截断的，延拓比较困难些。若只想要正确的求出各次谐波的幅值，这个有没有可行的方法呢？