请问王老师关于apFFT问题

a.gain

　　购买您的APFFT书，实际操作了下发现频偏校正和相位在对信号扫描时有点问题，如下程序：
　　clear all;clc;close all;
　　N=512;
　　t=-N+1:N-1;
　　f1=62.5;A1=0.87654;ph1=10.231;dc=0.143456;
　　f2=75.345;A2=0.0;ph2=60.567;
　　fs=1000;
　　FSX=fs/N;
　　y=A1*sin(2*pi*t*f1/fs+ph1*pi/180)+A2*sin(2*pi*t*f2/fs+ph2*pi/180)+dc;
　　y1=y(N:2*N-1);
　　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);
　　win2=winn/sum(winn);
　　y22=y2.*win2;
　　y222=y22(N:end)+[0 y22(1:N-1)];
　　y2_fft=fft(y222,N);
　　a2=abs(y2_fft);
　　p2=mod(phase(y2_fft)*180/pi,360);
　　ee=mod((p1-p2)/180/(1-1/N),1);
　　aa=(a1.^2)./a2*2;
　　r=round(f1/FSX);
　　fff=(floor(f1/FSX)+ee(r+1))*FSX
　　aaa=aa(r+1)
　　ppp=p2(r+1)
　　结果如下：fff =
　　64.4531
　　aaa =
　　0.8765
　　ppp =
　　280.2310
　　频率出错是校正多加了一根谱线位置，初相出错还没找到原因，望老师指教!谢谢!

a.gain

还有下面这个比值校正法求出相出入为什么很大，求其他两参数都非常精确比APFFT好
clear all;clc;close all;
N=1024;
t=0:2*N;
f1=50.1;a1=0.876541234;ph1=10.231;dc=0.143456;
f2=65.345;a2=0.0;ph2=60.567;
fs=1024;
FSX=fs/N;
y=a1*sin(2*pi*t*f1/fs+ph1*pi/180)+a2*sin(2*pi*t*f2/fs+ph2*pi/180)+dc;
y1=y(N:2*N-1);
y11=y1;
y11_fft=fft(y11,N);
ya=abs(y11_fft);
yp=mod(phase(y11_fft)*180/pi,360);
[max_a,k]=max(ya);
Rwk=real(y11_fft(k))/4-(real(y11_fft(k-1))+real(y11_fft(k+1)))/6+(real(y11_fft(k-2))+real(y11_fft(k+2)))/24;
Iwk=imag(y11_fft(k))/4-(imag(y11_fft(k-1))+imag(y11_fft(k+1)))/6+(imag(y11_fft(k-2))+imag(y11_fft(k+2)))/24;
kk=k+1;
Rwkk=real(y11_fft(kk))/4-(real(y11_fft(kk-1))+real(y11_fft(kk+1)))/6+(real(y11_fft(kk-2))+real(y11_fft(kk+2)))/24;
Iwkk=imag(y11_fft(kk))/4-(imag(y11_fft(kk-1))+imag(y11_fft(kk+1)))/6+(imag(y11_fft(kk-2))+imag(y11_fft(kk+2)))/24;
Zwk=sqrt(Rwk^2+Iwk^2);
Zwkk=sqrt(Rwkk^2+Iwkk^2);
aa=Zwk/Zwkk;
rr=(3-2*aa)/(1+aa);
f=(k-1+rr)*fs/N;
A=(2*pi*Zwk*rr*(1-rr*rr)*(4-rr*rr))/2/sin(rr*pi)/N*2;
mod(atan(Iwk/Rwk)*180/pi-rr*180, 360)

zhwang554

你将sin改为cos相位就对了, fft的相位谱对应cos的相位
你将24.5改为24.45频率就对了, 24.5正好对应校正数为1, 须加一个判断, 不是24.5, 就没有这个问题

[ 本帖最后由 zhwang554 于 2009-10-11 06:10 编辑 ]

a.gain

我看了校正是接近于1但实际不等于1，怎么判断呢？望指教

zhwang554

下面是程序是
whoru 发表于 2009-6-23 16:33   的贴子"求教：apFFT频谱校正问题！'中的程序
http://forum.vibunion.com/thread-83661-1-1.html
其中有判断语句,虽用於apfft/apfft法,也可用於fft/apfft法

function [fffa,AA,PH]=apFFT_Corret(k0,F2,F1,N);
clc;clear;clf;close all;
Fs=1000;
N=512;
t=-N+1:N*2-1;
f=62.50000;
A=0.87654;
Ph=10.231;
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 =62.50000000000000
校正的相位
PH =10.23099999999547
校正的振幅
AA =0.87654000000018

xiaokang

请教两个问题，1、如果我的信号事先不知道是什么样，具有几个频率值也未知，那么该程序中的CNum怎么取值（我随取了几个，就得到几个峰的值，很准）；2、如果我想画校正后的频谱图怎么画呢？（可否CNum取整个的点数，然后用plot（fff）？？？？我试了一下不对）？？？请教了，先谢谢！

xiaokang

我的目的就是对一个时间序列画频谱，点数已经确定，不能增加了，也不想补零（有泄漏），想得到更加准确的频谱图，主要就是所有的频率和幅值更准确（或者最大值的频率和幅值更准确，我在图中标出来），相位无所谓，采用什么频谱校正发好呢？谢谢

[ 本帖最后由 xiaokang 于 2009-10-27 06:05 编辑 ]

VibrationMaster

我觉得最简单的理论还是这个最简明
http://www.sciencedirect.com/sci ... 6adf18e13ad6ff8bc63

zhwang554

　　图贴不上
　　请参看在振动论坛下网址 日志<校正后的频谱图怎么画?>
　　http://forum.vibunion.com/UChome ... -blog-id-17964.html
　　[ 本帖最后由 zhwang554 于 2009-10-27 23:11 编辑 ]

xiaokang

王老师，在上面的网址上我没有权限发帖，只好又到这里，见谅！还想请教：那我最终想画横轴校正后的频率，纵轴校正后的幅值，这个实现不了吗？用plot(f,aa(1:N/2)),f是N转化来的，这样也不对啊，应该怎么画呢？谢谢

zhwang554

如要在f=1:N/2=1:256,频间为1赫的谱上标出校正后的频率为12.47赫,振幅为10的谱线, 将f 轴扩大100倍,f=1:25600,在x=1247处画振幅10

[ 本帖最后由 zhwang554 于 2009-10-29 23:19 编辑 ]

bocaoboluo

怎么计算apFFT的幅值呢？还是没有看明白，与FFT结果根号的关系吗？

zhwang554

　　由插值公式, apfft的泄漏是fft的泄漏的平方
　　h=2*pi*g.*(1-g.*g)./sin(pi*g)
　　其中g为频偏值
　　hanning窗的fft插值公式 S=A/h
　　所以 A=S*h (A为信号振幅值)
　　hanning窗的apfft插值公式 Sap=A/(h^2)
　　所以 A=Sap*h^2
　　5楼程序中除2是因插值公式对exp信号,对余弦信号幅值减半
　　[ 本帖最后由 zhwang554 于 2009-11-18 03:04 编辑 ]

bocaoboluo

我编写的一个函数为
x(i)=2*cos(2.2*pi*2*(i-127)/128+10*pi/180)，如果改成
x(i)=2*cos(2.2*pi*2*i/128+10*pi/180)  相位结果就不对了，为什么呢？
我自己编写算法，每个包含x(0) 的组合移位 ——求和——DFT的结果 与 每个包含x(0)的组合移位——每个组合DFT——每组FFT结果对应求和得出的结果不同，
问题：这种DFT结果对应频率点相加取平均 是不是等价于 每个数据先叠加再取平均呢？谢谢解答

bocaoboluo

我编写的一个函数为
x(i)=2*cos(2.2*pi*2*(i-127)/128+10*pi/180)，如果改成
x(i)=2*cos(2.2*pi*2*i/128+10*pi/180)  相位结果就不对了，为什么呢？
我自己编写算法，每个包含x(0) 的组合移位 ——求和——DFT的结果 与 每个包含x(0)的组合移位——每个组合DFT——每组FFT结果对应求和得出的结果不同，
问题：这种DFT结果对应频率点相加取平均 是不是等价于 每个数据先叠加再取平均呢？谢谢解答