请问王老师关于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
频率出错是校正多加了一根谱线位置，初相出错还没找到原因，望老师指教!谢谢!
还有下面这个比值校正法求出相出入为什么很大，求其他两参数都非常精确比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)
回复楼主a.gain的贴子
你将sin改为cos相位就对了, fft的相位谱对应cos的相位你将24.5改为24.45频率就对了, 24.5正好对应校正数为1, 须加一个判断, 不是24.5, 就没有这个问题
[ 本帖最后由 zhwang554 于 2009-10-11 06:10 编辑 ] 我看了校正是接近于1但实际不等于1，怎么判断呢？望指教
回复地板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
回复 5楼 zhwang554 的帖子
请教两个问题，1、如果我的信号事先不知道是什么样，具有几个频率值也未知，那么该程序中的CNum怎么取值（我随取了几个，就得到几个峰的值，很准）；2、如果我想画校正后的频谱图怎么画呢？（可否CNum取整个的点数，然后用plot（fff）？？？？我试了一下不对）？？？请教了，先谢谢！回复 5楼 zhwang554 的帖子
我的目的就是对一个时间序列画频谱，点数已经确定，不能增加了，也不想补零（有泄漏），想得到更加准确的频谱图，主要就是所有的频率和幅值更准确（或者最大值的频率和幅值更准确，我在图中标出来），相位无所谓，采用什么频谱校正发好呢？谢谢[ 本帖最后由 xiaokang 于 2009-10-27 06:05 编辑 ] 我觉得最简单的理论还是这个最简明
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
回复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 编辑 ]
回复 9楼 zhwang554 的帖子
王老师，在上面的网址上我没有权限发帖，只好又到这里，见谅！还想请教：那我最终想画横轴校正后的频率，纵轴校正后的幅值，这个实现不了吗？用plot(f,aa(1:N/2)),f是N转化来的，这样也不对啊，应该怎么画呢？谢谢回复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 编辑 ] 怎么计算apFFT的幅值呢？还是没有看明白，与FFT结果根号的关系吗？
回复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 编辑 ]
我编写的一个函数为
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结果对应频率点相加取平均 是不是等价于 每个数据先叠加再取平均呢？谢谢解答
回复 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