我要加入 登录
声振论坛 返回首页

bcyangbc的个人空间 http://home.vibunion.com/?85116 [收藏] [复制] [分享] [RSS]

日志

[转]全相位FFT算法的实现

已有 1285 次阅读2010-1-3 17:41 |

关键字:apFFT 频谱校正 全相位

传统FFT利用三角函数的正交性,将信号分离出来,从而将时域的信号变换到频域。但是,它有一个很重要的前提:输入的序列必须是周期内等间隔采样的值,这样,FFT计算的结果才是我们想要的。
实际的情况是,很难做到等间隔采样。比如,交流电的频率是变化的,并不是固定的50Hz。如果采用按照50Hz的信号来采样,则计算结果将无法反映原始信号。

为什么会出现上述的情况呢?因为FFT是将截断的序列做周期延拓而得到无限长的序列的。当不是等间隔采样时,做周期延拓后,在首尾相接的地方,就会出现信号的跳变,与原始信号不一致,自然也就得不到想要的结果,在频谱上的表现就是频谱出现泄露。

传统FFT的结果可以通过一些算法实现频谱校正,如全能重心法、比值法等等。但都是基于FFT的结果,精度有限。

全相位FFT算法是天津大学的王兆华和候正信教授提出的,具有初始相位不变和有效防止频谱泄露的特性,在他的相关书籍上做过严格的证明,并有相应的matlab程序。

我的项目中,需要分析50±0.5Hz工频信号,所以,研究apFFT算法,希望解决非同步采样的频谱校正问题。

apFFT的算法是先对数据进行全相位预处理,然后进行传统的FFT运算。FFT运算大家都熟悉,在这里,我说说全相预处理的过程。

全相位预处理过程如下:
1,构成一个N点的汉宁窗。
2,汉宁窗对自己求卷积,得到2N-1点的卷积窗。
3,求2N-1点的卷积窗的和。
4,将卷积窗的每一项除以卷积窗的和,得到2N-1点的归一化卷积窗。
5,将数据的1:2N-1项和归一化卷积窗相乘,得到加窗的2N-1项。
6,将第1项和N+1项,第2项和N+2项 ... 第N-1项和第2N-1项相加,得到经过全相预处理的N点序列。
接下来就只需要进行FFT,就可以得到apFFT的结果。

如果我们观察结果全相位预处理的信号,由原来的带有跳变的信号变为平滑连续的信号,使得周期延拓后不会出现信号的跳变。

下面是apFFT的matlab程序:
close all;clc;clear all;
N=256;
t=-N+1:2*N-1;
f1=39.8;
f2=51.345;
A1=5;
ph1=30;
s=A1*cos(1*(2*pi*t*f1/N/50-ph1*pi/180));%+A1*cos(1*(2*pi*t*f2/N+ph1*pi/180));
win=hanning(N)';
win1=hanning(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);
y1=s1.*win2;
y1a=y1(N:end) + [0 y1(1:N-1)];
Out1=fft(y1a,N);
a1=abs(Out1);
p1=mod(phase(Out1),2*pi);
s2=s(1+N:3*N-1);
y2=s2.*win2;
y2a=y2(N:end) + [0 y2(1:N-1)];
Out2=fft(y2a,N);
a2=abs(Out2);
p2=mod(phase(Out2),2*pi);
g=mod((p2-p1)/pi/2,1);
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a2)/2;
rr=round(f1);
disp('频率校正值')
fff=floor(f1)+g(rr+1)
disp('振幅校正值')
aaa=aa1(floor(f1)+1);
disp('初相位校正值')
ppp=p1(rr+1)*180/pi;

通过在VC下编码验证,apFFT具有很高精度的初始相位。采用apFFT/apFFT时移相位差法,可以高精度的修正频率偏差和幅值偏差。这个算法对于非同步采样处理的意义,不言而喻。与传统fft频谱校正的主要区别在于,apFFT的初始相位是准确的,相位差也就是准确的,频率也就可以校正准确,所以具有更高的精度。王教授对各种校正方法做过详细的测试和对比,若希望了解相关内容,可查阅《数字信号全相位谱分析与滤波技术》一书

发表评论 评论 (1 个评论)

回复 李元 2010-1-5 13:33
向大侠学习!

facelist doodle 涂鸦板

您需要登录后才可以评论 登录 | 我要加入

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-5 23:50 , Processed in 0.122088 second(s), 16 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部