a.gain 发表于 2009-10-10 23:37

请问王老师关于apFFT问题

本帖最后由 wdhd 于 2016-9-20 13:05 编辑

  购买您的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)+;
  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 发表于 2009-10-10 23:44

还有下面这个比值校正法求出相出入为什么很大,求其他两参数都非常精确比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(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 发表于 2009-10-11 06:06

回复楼主a.gain的贴子

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

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

a.gain 发表于 2009-10-11 09:10

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

zhwang554 发表于 2009-10-11 09:39

回复地板a.gain的贴子

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

function =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)+;
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 =62.50000000000000
校正的相位
PH =10.23099999999547
校正的振幅
AA =0.87654000000018

xiaokang 发表于 2009-10-27 05:44

回复 5楼 zhwang554 的帖子

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

xiaokang 发表于 2009-10-27 06:01

回复 5楼 zhwang554 的帖子

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

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

VibrationMaster 发表于 2009-10-27 08:50

我觉得最简单的理论还是这个最简明
http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6TJ5-4V59VYT-4&_user=1479053&_coverDate=06%2F30%2F2009&_alid=1065213591&_rdoc=2&_fmt=high&_orig=search&_cdi=5301&_sort=d&_docanchor=&view=c&_ct=19&_acct=C000053031&_version=1&_urlVersion=0&_userid=1479053&md5=edc3641aa56ec6adf18e13ad6ff8bc63

zhwang554 发表于 2009-10-27 23:10

回复6楼Xiaokang 的贴子----校正后的频谱图怎么画?

本帖最后由 wdhd 于 2016-9-20 13:07 编辑

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

xiaokang 发表于 2009-10-29 14:44

回复 9楼 zhwang554 的帖子

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

zhwang554 发表于 2009-10-29 20:49

回复10楼Xiaokang 的贴子

如要在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 发表于 2009-11-17 16:26

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

zhwang554 发表于 2009-11-17 23:10

回复12楼bocaoboluo 的贴子

本帖最后由 wdhd 于 2016-9-20 13:07 编辑

  由插值公式, 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 发表于 2009-11-18 22:34

我编写的一个函数为
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 发表于 2009-11-18 22:34

回复 13楼 zhwang554 的帖子

我编写的一个函数为
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结果对应频率点相加取平均 是不是等价于 每个数据先叠加再取平均呢?谢谢解答
页: [1] 2
查看完整版本: 请问王老师关于apFFT问题