快速傅里叶变化

快速傅里叶变化（共6篇）

快速傅里叶变化 篇1

1 音乐频谱显示介绍

音乐频谱显示, 就是通过硬件或软件的方法, 随着音乐强弱和节奏的变化, 将音频信号中每个频率成分的幅值信息通过屏幕实时地显示出来。

2 硬件设计

其大体设计思路为:首先对音频信号进行采集, 经过数据处理后输出至端口, 最后驱动屏幕显示。最小系统, 是指芯片能够正常工作所需的最小电路。在最小系统的基础上进行外围电路的扩展, 便构成了整个系统电路。它外围电路主要包括电源电路、端口扩展电路和显示电路。

系统采用12V交流电输入, 经过整流、滤波、降压后, 将电源分为+9V、+5V、+3.3V三个部分, +9V电源主要为+3.3V电源模块提供输入电压, +3.3V主要为lm3s615提供电源, +5V为显示屏提供电源。

lm3s615有34个IO口, 在本系统中, 屏幕为18*20的点阵屏, 其中20根管脚用于传输数据, 18根管脚用于每列的控制。数据端采用串行输入并行输出芯片74HC595, 它最少只需3个IO口即可实现多位数据的输出, 而总的时延仅为几十纳秒。控制端采用3-8译码器74LS138。为满足本系统的要求, 使用了三块74HC595和三块74LS138来对端口进行扩展。最终使用了7个IO口, 即可实现对38个管脚的控制, 极大地节省了IO口的使用。

显示将使用20个三极管, 来控制LED的亮灭。LED点阵使用了360个普通的发光二极管, 按行20个和列18个依次排布。

3 软件设计

单片机以41 KHz的速率采集到的音频信号都是时域里面的离散信号, 需要对数据进行离散傅里叶变换 (DFT) 。如果直接用DFT算法进行谱分析和信号的实时处理, 那是不切实际的。考虑到单片机运算速度有限的特点, 在本系统中采用了分裂基快速算法, 使运算效率进一步提高, 主要是对蝶形运算和倒序运算的编程。为保持频谱显示的协调性, 加入了峰值缓慢下降程序。

3.1 蝶形运算

完成一个蝶形运算, 需要一次复数乘法和两次复数加法。经过一次分解后, 计算1个N点DFT共需要计算两个N/2点DFT和N/2个蝶形运算。而计算一个N/2点DFT需要 (N/2) 的平方次复数乘法和N/2 (N/2-1) 次复数加法。由此可见, 仅仅经过一次分解, 就使运算量减少近一半, 所以可以对N/2点DFT作进一步分解。从而使运算量大大地减少。可以归纳出一些对编程有用的运算规律:第L级中, 每个蝶形的两个输入数据相距B=2L-1个点;每级有B个不同的旋转因子;同一旋转因子对应着间隔为2L点的2M-L个蝶形。

总结上述运算规律, 便可采用下述运算方法。先从输入端 (第1级) 开始, 逐级进行, 共进行M级运算。在进行第L级运算时, 依次求出B个不同的旋转因子, 每求出一个旋转因子, 就计算完它对应的所有2M-L个蝶形。这样, 可用三重循环程序实现DIT-FFT运算。根据基2-FFT算法, N点FFT运算可以分成log2N级, 每一级都有N/2个蝶形运算。可见每个蝶形运算的输出都是由其输入值与某一正弦函数和余弦函数的乘积累加得到的。

3.2 倒序运算

以8点为例, 来阐述倒序它的规律所在。在运算M级蝶形运算之前应先对序列进行倒序。为叙述方便, 用J表示当前倒序数的十进制数值。对于N=2M, M位二进制数最高位的十进制权值为N/2, 且从左向右二进制位的权值依次为N/4, N/8, …, 2, 1。

因此, 最高位加1相当于十进制运算J+N/2。如果最高位是0 (J<N/2) , 则直接由J+N/2得下一个倒序值;如果最高位是1 (J≥N/2) , 则先将最高位变成0 (J⇐J-N/2) , 然后次高位加1 (J+N/4) 。但次高位加1时, 同样要判断0、1值, 如果为0 (J<N/4) , 则直接加1 (J⇐J+N/4) , 否则将次高位变成0 (J⇐J-N/4) , 再判断下一位;依此类推, 直到完成最高位加1, 逢2向右进位的运算。

3.3 峰值缓慢下降

为避免数据改变过于频繁, 使显示故加此算法其主要思路是:当前采样值与前一次采样值相比较, 若比前一次值高, 则更新数据, 否则, 在前一次值的基础上缓慢下降, 直至为零。中途若有数据更高, 则再次更新数据。如此循环, 便可实现峰值的缓慢下降。

4 系统测试

在完成硬件和软件的设计后, 需对整个系统进行测试以保证方案的可行性。测试可分为FFT算法的测试、数据采集测试和整体测试。

FFT算法测试:该测试为验证FFT算法得出的结果。主要测试环境为Microsoft VisualC++6.0。

数据采集测试:将音频信号输入至单片机, 再通过串口观察数据的变化, 观察工具为windows xp的超级终端。

整体测试:该测试主要是观察输出的频谱是否与电脑上显示的频谱一致。第一次输入1 KHz的正弦信号, 第二次输入普通音频信号, 再分别作对比。通过比较, 得到的结果与预期结果基本一致, 系统测试通过, 设计基本完成。将该频谱显示器放置在音响旁, 即实惠, 又能为我们的音乐增添一点色彩。

参考文献

[1]高西全, 丁玉美.数字信号处理[M].西安:西安电子科技大学出版社, 2008 (08) :11-223.

[2]杨颂华, 冯毛官, 孙万蓉, 胡力山.数字电子技术基础[M].西安:西安电子科技大学出版社, 2000 (07) :2-45.

[3]张富.C及C++程序设计[M].北京:人民邮电出版社, 2005 (07) 67-189.

[4]籍顺心.单片机的C语言应用程序设计 (第四版) [M].北京:北京航空航天大学出版社, 1999 (09) :2-129.

快速傅里叶变化 篇2

关键词：快速傅里叶变换,逆序,快速算法

0 引言

快速傅里叶变换(FFT)是数字信号处理中常用的方法,在FFT算法中都要涉及到逆序数的计算,本文称为整序。目前用于FFT中整序的算法较多,较为常见的为大多数文献上列出的“逢二退一”法[1,2,3],该方法编程简单,内存占用少,但运算量不是最优;文献[4]对“逢二退一”法进行了优化,运算量有所减小;文献[5,6]提出的生成法,较 “逢二退一”法有效的减少了运算量。本文介绍一种新方法,较生成法更进一步的减少了运算量,提升了运算速度,而且编程也很简单。

1 算法原理

首先,本文约定二进制数从左到右为从高位到低位。表1是点数N=16的原序数(十进制)及其逆序数(二进制)的对照关系。

观察表1可知:0的逆序数一定为0000,1的逆序数只要将0的逆序数最高位由0变为1(即一次加法运算)得到,继续观察可知,原序数2,3的逆序数分别可由原序数0,1的逆序数在次高位由0变为1得到,进一步观察,原序数4,5,6,7的逆序数分别可由原序数0,1,2,3的逆序数在第三高位由0变为1得到,原序数8至15的逆序数分别可由原序数0至7的逆序数在第四高位由0变为1得到。

以上观察结果可以推广到N=2M的情况,其规律可用以下定理描述:

定理:设N=2M,i=1,2,…,M,Yj为原序数j的逆序数,Y0=0,在N个逆序数中,原序数为2i-1, 2i-1+1, 2i-1+2,…,2i-1的逆序数可分别由原序数为0,1,2,…,2i-1-1的逆序数加2M-i得到,即

Y(2i+1+j)=Yi+2M-I,j=0,1,2,……,2i-1-1

2 算法分析

本文介绍的整序快速实现算法中,在欲求的N(N=2M)个逆序数中,前2i(1<=i<=M)个逆序数的后半部分可由前半部分逆序数依次通过一次加法生成,此过程共需要2i-1次加法,这样通过M次迭代(i从1到M),只需要N-1次加法运算即可生成N个逆序数,即

${\rm N}-1=2{}^{{\rm M}}-1={\displaystyle \sum _{i=0}^{{\rm M}-1}2}{}^i$

为了便于比较,下面简单介绍生成法[5,6]的基本原理:

在生成法中,只要知道2 的i-1阶逆序数,即可求出2 的i阶逆序数。2 的i阶逆序数的前半部分是由i-1阶全部逆序数顺次乘2 组成,后半部分是前半部分顺次加1生成。此过程经过M次迭代(i从1到M)即可得到N(N=2M)个逆序数。

生成法算法中,2的i阶逆序数的前半部分分别需要一次乘2运算,即一次迭代共需要2i-1次乘2,M此迭代共需要N-1次乘2,即

${\rm N}-1=2{}^{{\rm M}}-1={\displaystyle \sum _{i=0}^{{\rm M}-1}2}{}^i$

后半部分中一个逆序数需要一次加法运算,即一次迭代共需要2i-1次加法,M此迭代共需要N-1次加法,即

${\rm N}-1=2{}^{{\rm M}}-1={\displaystyle \sum _{i=0}^{{\rm M}-1}2}{}^i$

所以,生成法一共需要N-1次乘2运算和N-1次加法运算。

几种方法的运算量比较如表2所示。

比较可知,本文提出的算法在运算次数上均优于其他两种算法。

本文中提出的算法在内存使用方面与生成法相同,当样点数为N时,需占用N单位内存(根据N值大小,单位可为字节、双字节、四字节),较逢二退一法使用较多的内存单元。

3 算法实现

根据以上公式,取Y0=0(0的逆序数一定为0),当i分别取1,2,….M时,即可循环得到2M个逆序数。

C语言实现如下:

运行后,Y[j]中存放原序数j对应的逆序数。

4 结束语

本文提出的算法使用的运算量较少,均小于参考文献中的算法的运算量。应用该算法会大幅度提高整序速度,进而提高整个FFT运算的速度。

参考文献

[1]程培青.数字信号处理[M].清华大学出版社,2001.

[2]丁玉美,高西金.数字信号处理[M].2版.西安电子科技大学出版社,2001.

[3]A.V奥本海姆,R.W谢弗.离散时间信号处理[M].黄建国,刘树棠,译.北京科学出版社,1998.

[4]高丽,刘卫新,张学智.FFT标准整序算法的优化[J].探测与控制学报,2004,26(2):62-64.

[5]张学智,沈虹.实现快速傅里叶变换中逆序的新方法[J].西安工业学院学报.

快速傅里叶变化 篇3

传统的傅里叶变换(FFT)对平稳信号的处理效果很好,但当信号频率随时间变化时,FFT就显得有些力不从心了。分数阶傅里叶变换(Fractional Fourier Transform,FRFT)可以很好地弥补FFT的不足,特别是处理线性调频信号(LFM)时,能够得到令人满意的结果。

FRFT也称为角度傅里叶变换(AFT) 或者旋转傅里叶变换(RFT),其定义式为:

$X{}_p(u)=R{}^{\alpha }[x(t)]={\displaystyle {\int }_{-\infty }^{\infty }{\rm K}}{}_p(t,u)x(t)\text{d}t(1)$

式中变换核取作:

${\rm K}{}_p(t,u)=\{\begin{array}{l}\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}(\frac{1}{2}u{}^2\cot \alpha -ut\text{c}\text{s}\text{c}\alpha +\frac{1}{2}t{}^2\cot \alpha )},\alpha \ne n\text{π}\\ \delta (t-u),\alpha =2n\text{π}\\ \delta (t+u),\alpha =(2n+1)\text{π}\end{array}$

则:

$\begin{array}{l}X{}_p(u)=\{F{}^{D}x\}(u)={\displaystyle {\int }_{-\infty }^{\infty }x}(t)k{}_p(t,u)\text{d}t\\ =\{\begin{array}{l}\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }{\displaystyle {\int }_{-\infty }^{\infty }x}(t)\text{e}{}^{\text{j}\frac{t{}^2}{2}\cos \alpha }\text{e}{}^{-\text{j}ut\text{c}\text{s}\text{c}\alpha }\text{d}t,\alpha \ne n\text{π}\\ x(t),\alpha =2n\text{π}\\ x(-t),\alpha =(2n+1)\text{π}\end{array}(2)\end{array}$

其中n为整数,即n∈Z。α=pπ/2称为分数阶Fourier变换的阶数,并有$R{}^{\alpha }=R{}^{p\frac{\pi }{2}}=F{}^p$。Kp(t,u)称为FRFT的核函数。Xp(u)称为x(t)的p阶Fourier变换。FRFT是一种线性算子,记为Fp,他满足以下性质:

(1) FRFT变换为线性算子;

(2) F0[x(t)]=F4[x(t)]=x(t)(恒等变换);

(3) F1[x(t)]=F5[x(t)]=X(ω)(标准Fourier变换);

(4) 广义Fourier变换算子为加性算子,即有Fp+q=FpFq。

2 采用分解方法计算FRFT的步骤

分数阶Fourier变换可以具体分解为以下三个主要的计算步骤:线性调频信号乘法;线性调频信号卷积;另一个线性调频信号乘法。假定p∈[-1,1],则我们可以将经过量纲归一化的信号f(x)的分数阶Fourier变换式(2)分解为以下三步运算:

$f{}_p(x)=\text{e}{}^{-\text{j}\text{π}x{}^2\tan (\alpha /2)}{g}^{\prime }(x)(3)$

和:

$\begin{array}{l}{g}^{\prime }(x)=A{}_{\alpha }{\displaystyle {\int }_{-\infty }^{\infty }\text{e}}{}^{\text{j}\text{π}\beta (x-x^{\prime }){}^2}g(x^{\prime })\text{d}{x}^{\prime }(4)\\ g(x)=\text{e}{}^{-\text{j}\text{π}x{}^2\tan (\alpha /2)}f(x)(5)\end{array}$

式中,α=pπ/2,β=csc α,而g(x)和g′(x)表示中间结果,并且:

$A{}_{\alpha }=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}(6)$

即是说,分数阶Fourier变换的数值计算的顺序如下:先计算式(式(5)),再计算式(4),最后计算式(3)。下面是每一步计算的有关细节。

第一步:将函数f(x)与线性调频函数相乘(式(5))。注意,g(x)的频率带宽与时间带宽乘积可以是f(x)的相应带宽乘积的两倍,所以要求g(x)的采样间隔为1/(2Δx)。如果f(x)样本值的采样间隔是1/Δx,那么就需要对这些样本值进行插值,然后再与线性调频函数的离散采样值相乘,以得到所希望的g(x)的采样。

第二步:将g(x)与一线性调频函数作卷积式(式(4))。注意,由于g(x)是带限信号,所以线性调频函数也可以用其带限形式代替而不会有任何影响。也就是说,我们可以取:

$\begin{array}{l}{g}^{\prime }(x)=A{}_{\alpha }{\displaystyle {\int }_{-\infty }^{\infty }\text{e}}{}^{\text{j}\text{π}\beta (x-x^{\prime }){}^2}g(x^{\prime })\text{d}{x}^{\prime }\\ =A{}_{\alpha }{\displaystyle {\int }_{-\infty }^{\infty }h}(x-x^{\prime })g(x^{\prime })\text{d}{x}^{\prime }(7)\end{array}$

式中:

$h(x)={\displaystyle {\int }_{-\text{Δ}x}^{\text{Δ}x}{\rm H}}(\gamma )\text{e}{}^{\text{j}2\text{π}\gamma x}\text{d}\gamma (8)$

这里:

${\rm H}(\gamma )=\frac{1}{\sqrt{\beta }}\text{e}{}^{\text{j}\text{π}/4}\text{e}{}^{-\text{j}\text{π}\gamma {}^2/\beta }(9)$

是函数ejπβx2的Fourier变换。于是,式(7)的离散形式为:

$g^{\prime }(\frac{m}{2\Delta x})={\displaystyle \sum _{n=-{\rm N}}^{{\rm N}}h}\left(\frac{m-n}{2\text{Δ}x}\right)g\left(\frac{n}{2\text{Δ}x}\right)(10)$

这一离散卷积可以利用快速Fourier变换计算。

第三步:计算式(3)得到f(x)的分数阶Fourier变换fp(x)的采样值$f{}_p\left(\frac{m}{2\text{Δ}x}\right)$。由于假定f(x)的所有变换都是带限的,他们位于区间$\left[-\frac{1}{2}\text{Δ}x,\frac{1}{2}\text{Δ}x\right]$,所以需要用因子2对$f{}_p\left(\frac{m}{2\text{Δ}x}\right)$进行二抽一采样,以得到离散采样$f{}_p\left(\frac{m}{2\text{Δ}x}\right)$。因此,对于不是π/2整数倍的角度,分数阶Fourier变换的计算对应以下步骤:

(1) 原信号与一线性调频函数相乘;

(2) Fourier变换(其变元乘以尺度系数csc α);

(3) 再与一线性调频函数相乘;

(4) 乘以一复幅值因子。

3 对信号的FRFT处理及仿真图

首先要给出一个输入信号x,然后根据分解方法编出FRFT快速算法,根据不同的p值输入信号x会生成不同的曲线,其中p∈(0,4)。找出每个p值时FRFT的输出最大值点,组成一个一维数组m1,再从m1中找出一个最大值,该值所对应的p值就是FRFT变换的最佳角度α=pπ/2。

例1 假设输入信号为方波C=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ],长度N=73,在p=0～1之间进行FRFT变换,仿真结果如图1所示。

由图1(d)中可以看出当p=1时,α=π/2就是普通的傅里叶变换,这也验证了分数阶傅里叶变换的正确性。

例2 输入信号为多项式的线性调频信号fmpoly(N,p),x=exp(j*pi*(n-center)·^p/(N/2)^(p-1)/p),其中N为信号长度,p是信号的多项式幂,当p=2时该信号为chirp信号。取x=fmpoly(63,2)。

由图2可以看出,在图2(e)中p=1.5的时候,即α=3π/4时形成了一个冲击信号,说明了在此角度上信号的能量最好地集聚在一点上,由此可以识别出信号的调频系数,检测出信号的参数,这就是FRFT处理LFM信号的显著作用。

4 结 语

分数阶傅里叶变换是近二十年来发展起来的一种全新的信号时频分析工具,在很多方面得到了十分广泛的应用。而其快速算法的研究则对扩展其应用领域有着十分重要的意义。本文提出了一种有效并能准确计算FRFT 的新算法。该算法具有易实现、易理解、精度较高等优点,相信FRFT将会受到更广泛的重视,在信号处理领域会有良好的应用前景。

参考文献

[1]平先军,陶然,周思永,等.一种新的分数阶傅里叶变换快速算法[J].电子学报,2001,29(3):406-408.

[2]Lufs B Almeida.The Fractonal Fourier Transform andTime-Frequency Representations[J].IEEE Transactionson Signal Processing,1994,42(11).

快速傅里叶变化 篇4

关键词：心律失常,快速傅里叶变换,归一化相关系数,模板匹配

0 引言

心血管疾病严重威胁着人类健康。心血管疾病, 尤其是心律失常的检测依赖于心电图。实现心电图 (electrocardiogram, ECG) 的智能分析和自动识别, 对心血管疾病的辅助诊断具有重要的临床意义。ECG智能分析技术包括心电信号预处理、波形检测、心电诊断等。预处理重在对含有噪声的心电信号消除干扰、滤除噪声;波形检测的目的在于寻找信号特征点, 提取特征参数;而心电诊断技术则主要依据心电波形特征点的识别结果, 对信号进行分类。其中ECG的自动识别方法包括特征提取、模板匹配、支持向量机、人工神经网络等[1,2,3,4]。

模板匹配算法在ECG自动识别中具有较高的应用价值, 心电模板匹配算法包括线性预测算法和波形形态分类算法等。线性预测算法采用的基本思想是通过多组ECG样本数据的线性组合来估计某一未知ECG数据;波形形态分类算法则比较待测心电与原规定 (或学习到) 的心电模板或规则的近似程度, 据此判断待测信号与模板是否属于同一类别。基于信号波形形态特征分析的算法, 能对具有相似特征参数 (如QRS复波幅度、宽度等) 不同形态特征的心电进行区分, 与医生进行人工诊断时所使用的方法具有一致的思想[5,6,7]。本文采用基于波形形态分类的模板匹配算法来进行心律失常研究。

1 理论基础

1.1 归一化相关系数

ECG信号与模板信号的相似程度可通过计算它们在指定区间的相关性来判定[8]。本文研究利用归一化相关系数来衡量待测信号与ECG模板的相似程度。

设有2个离散信号x (n) 和y (n) , 将其中一个信号扩展a倍, 则二者之间的误差能量ε2为:

为了比较2个信号的相似程度, 须使误差能量最小, 故有:

从而得到a以及最小误差能量:

以x (n) 的能量为基准, 得到相对最小误差能量e2min为:

令

则ρxy称为归一化相关系数或相干系数。在一个序列移动的情况下, 相干系数变成相干函数, 以ρxy (m) 表示:

又由离散信号x (n) 和y (n) 的相关函数:

可得:

相干函数ρxy (m) 即为归一化至[-1, 1]范围内的相关函数[9], 若m=0, 则ρxy为通常意义下的互相干系数。由于相关函数rxy (m) 的取值与x (n) 和y (n) 的序列大小有关, 难以比较各组信号之间的相关程度;而相干函数值的大小与信号值的大小无关, 不仅可以表示同相相干性, 还可以表示移相 (或延时) 相干性, 因此, 选取相干函数作为模板匹配的判断依据。

1.2 利用快速傅里叶变换 (fast Fourier transform, FFT) 计算归一化相关系数

2 个长为N的离散时间序列x (n) 和y (n) 的互相关函数为:

其离散傅里叶变换为:

其中, X (k) =DFT[x (n) ], Y (k) =DFT[y (n) ], Rxy (k) =DFT[rxy (τ) ]

将x (n) 、y (n) 的离散傅里叶反变换代入互相关函数公式, 得:

因x (n) 是实序列, 所以x (n) =x* (n) , 得

由

得

当x (n) =y (n) 时, 得x (n) 的自相关函数为:

由此可以看出, 互相关和自相关函数的计算可通过FFT变换求得, 而相干函数和相干系数均与相关函数有关, 故可通过FFT变换实现归一化相关函数和归一化相关系数的计算。

2 研究方法

2.1 数据来源

本研究所采用的心电数据来自MIT-BIH心律失常数据库, 数据库中的每条记录由3个文件组成:头文件 (hea) 、数据文件 (dat) 和注释文件 (atr) 。头文件包括对数据存储格式的说明、信号的通道数、导联方式、采样精度和长度、检验位以及患者的个人信息记录, 数据文件为该患者的心电数据, 注释文件标注了每例信号所记录的时间段以及类型[10]。待测心电为从MIT-BIH心律失常数据库中挑选的正常或异常心拍, 模板心电则采用同一采集者的同一导联的正常心拍。文中采用编号为202和213患者的多个心拍记录数据, 模板信号和待匹配的心电拥有相同长度, 均为一个心动周期。

2.2 算法设计

传统相关函数的计算采用直接计算卷积的方法实现, 但卷积的运算量大、耗时多、实时性差。由上述理论推导可知, 互相关和自相关函数的计算可通过FFT求得, 因此本文先利用计算相关函数, 再利用相关函数来求相干函数和相干系数, 进而得到互相干系数。此方法可以大幅减少计算量, 提高处理速度。具体算法流程如图1所示。

首先利用FFT计算待测信号和模板信号的自功率谱和互功率谱, 然后再进行FFT反变换, 求得待测信号和模板信号的自相关函数和互相关函数, 最后利用相关函数得到互相干函数, 进而得到归一化相关系数。将计算得到的待测信号与模板信号的互相关系数与预设阈值比较 (本文使用经验值0.92) , 若此值大于预设阈值, 表明2个信号匹配成功, 判定该待测ECG信号与匹配模板是同一类别, 为正常信号;否则, 判为异常。

3 结果与讨论

选取数据库中记录为202和213的多个心拍数据进行仿真实验, 其中通道1为肢体Ⅱ导联, 通道2为胸导联V1。

(1) 模板信号为正常心拍, 待测信号为同一个体的异常心拍, 得到结果见表1。

由表1可以看出, 以正常心拍作模板, 得到2通道 (V1导联) 异常心拍的判断正确率为100%, 1通道有2例心室融合心跳的异常心拍被误判为正常。

(2) 模板信号为正常心拍, 待测信号为同一个体的正常心拍, 交叉组合选取正常模板与待测信号共12组, 对待测心拍进行判别, 结果见表2。

由表2可以看出, 1通道 (肢体Ⅱ导联) 的正常心拍判断准确率为100%, 胸导联V1中有3例正常心拍被误判。

由以上2组对比实验可知, 以正常心拍作模板, 运用基于FFT的心电模板匹配算法, 对肢体导联Ⅱ导联的正常信号能够全部检出, 异常信号有少部分被错判为正常;对胸导联V1的异常信号能够全部检出, 正常信号有一半被误判。因此, 对于心律失常的检测, 建议以胸导联V1信号的识别结果为依据;而对于健康体检, 则建议以肢体导联Ⅱ导联的结果为准;或者同时采集V1导联和Ⅱ导联的信号进行2次检测, 综合判断结果, 提高检测的准确率。

本文所选取的信号来源于MIT-BIH数据库, 如果能够从医院临床取得心电数据进行研究, 可能能够得到更加符合我国国情的统计结果, 下一步拟从周边医院开展调研、取得数据并进行深入研究。

4结语

利用FFT算法来计算归一化相关函数, 有可靠的理论基础;通过计算归一化相关系数来进行心电模板匹配的算法简便、快捷、准确率较高, 在心律失常疾病的辅助诊断中有一定的临床价值。但此算法要求待测信号和模板信号拥有相同的数据长度, 所以信号的采集受到限制, 不利于来源不同的各种心电信号之间的匹配。因此, 寻求能够实现可变模板自动识别的模板匹配算法将会有更加重大的意义。

参考文献

[1]周群一, 吕旭东, 段会龙.ECG心搏模式识别[J].生物医学工程学杂志, 2005, 22 (1) :202-206.

[2]杨煦, 张跃.基于支持向量机的心律失常识别方法[J].计算机工程与设计, 2007, 28 (18) :4 442-4 445.

[3]刘春玲, 王旭.基于改进的自适应小波神经网络的心电信号分类[J].系统仿真学报, 2007, 19 (14) :3 281-3 282.

[4]王鸿鹏, 戴博, 杨孝宗.心电自动分析系统[J].北京生物医学工程, 2007, 26 (4) :387-389.

[5]Hsia P W, Jenkins J M, Shimoni Y, et al.An automated system for ST-segment and arrhythmia analysis in exercise radionuclide ventriculography[J].IEEE Trans on BME, 1986, 33 (6) :585-592.

[6]程小明, 林金森, 张正国.高分辨心电图中模板匹配算法的改进[J].中国生物医学工程学报, 1999, 18 (1) :89-96.

[7]Lin K P, Chang W H.QRS feature extraction using linear prediction[J].IEEE Trans on BME, 1989, 36 (10) :1 050-1 055.

[8]张集墨, 张跃, 杨波.室性早搏的多模板匹配自适应识别算法[J].计算机工程与设计, 2011, 32 (8) :2 885-2 888.

[9]饶妮妮.生物医学信号处理[M].成都:电子科技大学出版社, 2006:45-46.

快速傅里叶变化 篇5

1 快速傅里叶变换和小波变换原理

1.1 快速傅里叶变换

快速傅里叶变换原理:在DFT中, 我们令系数WN=e-j2π/N, 由此可看出系数WN的一些性质。

为简单起见, 我们取N为2的整数次幂, 根据系数的对称性, 则有:

公式 (1) 只计算出了前N/2个点, 根据系数的对称性, 则有:

此时, 可以将1个N点的变换分解为2个N/2点的变换, 并且可以依据这种模式继续分解下去。这就是Cooley-Turkey的FFT算法的基本原理。它基本上分为时间抽取 (DIT-FFT) 算法和频率抽取 (DIF-FFT) 算法两类。

1.2 小波变换

“小波变换”概念最早于1984年由J.Morlet提出, 其基本思想是把信号投影在由一簇基函数张成的空间上。利用小波分析, 不仅能将信号在时间和频率上独立分解, 还能保证不丢失原有的信号特征, 被誉为信号分析中的“显微镜”。

2 风机振动的监测和分析

在整个风机系统中, 电机和风机轴承是核心部件, 也是风机运行故障的主要来源, 因此, 应被当作监测对象。据统计, 在所有风机故障中, 近70%的故障与转轴及其组件系统有直接的关系。因此, 对于风机而言, 检测点最好设在轴承部位, 且选择探头与机械接触较好、刚度较高之处作为测试点。

2.1 监测参数的选择

监测参数的选择原则为:对于低频 (振动频率小于10 Hz) 振动, 常取位移作为测量参数;对于中频 (振动频率在10~1 000 Hz之间) 振动, 取速度作为测量参数;对于高频 (振动频率在1~10 k Hz) 振动、随机振动等, 常将加速度作为测量参数。本文选择振动速度作为测量参数。

2.2 数据采集

振动一般由一系列简谐振动分量、其他形式的振动及噪声叠加形成, 因而对于振动信号的监测, 通常选取振动加速度、振动速度或振动幅值等作为测量参数。依据振动参数选择原则, 本文选用SG-2磁电式速度传感器。该速度传感器可输出微弱的电荷信号, 经电荷放大器和电压放大器后送入A/D转换器;将采集到的振动数据输入工具软件中进行信号数据分析和处理, 以此获取风机运行的振动状态及可能出现的故障点。

2.3 风机振动的判别标准

本文选用的是ISO 02372振动标准。ISO 02372振动标准是一种根据轴承振动烈度来评定机器质量的标准。本文所研究的风机属于第Ⅲ类风机。通常, 这类风机A振动区域的速度值为0.28~1.8 mm/s, B振动区域的速度值为1.8~4.5 mm/s, C振动区域的速度值为4.5~11.2 mm/s, D振动区域的速度值大于11.2 mm/s。根据此判别标准, 可确定设备的维修情况, 加强对C区域振动的监测。必要时, 还要加大监测力度。

3 基于FFT与小波变换的故障信号提取

3.1 风机技术参量

本文研究的风机技术参量:主风机安装轴承处的轴直径190mm, 叶片数13, 双列向心短圆柱轴承转速1 000 r/min, 转动频率fr=n/60=16.67 Hz, 叶轮通过频率fz=frZ=16.67×13=216.71 Hz (Z为叶片数) ;电动机转速1 000 r/min, 转频16.67 Hz。考虑到高次谐波, 转频的频段出现在中高频段, 选取振动速度作为测量参数。

3.2 故障信号的时域分析

在不平衡故障中, 主要振动特征就是存在以工频为主的重复性成分, 因而其时间波形表现出显而易见的正弦波形状, 振动信号表现为明显的正弦变化, 初步判断该设备的故障类型为装置不平衡故障。

3.3 故障信号的FFT分析

经过FFT分析获得的频谱分析图仅表现出某一个频率在所有采样振动信号中的总强度, 不能很好地反映该频率所对应的时间方面的具体信息, 即频谱分析能反映频率特征, 对频率的时间分辨率较低。

3.4 故障信号的小波分析

小波分析法是一种分辨率较高的时频分析法。运用这一方法不仅可以进行时域上的分析, 还可以进行频域上的分析。另外, 运用小波分析法不仅能精确定位短时高频信号, 还能准确分析低频信号。本文选Biorthogonal双正交样条小波作为小波基函数, 重构滤波器阶数Nr=6, 分解滤波器阶数Nd=8, 分解层数为7, 提高了频率范围域中的分辨率和分析时间域中的精度。

3.5 故障诊断

引起风机异常的原因为转子组件不平衡。考虑到风机恶劣的工作环境及工作过程中的突发因素, 初步判断是因风机在运转过程中, 其转子出现一定的磨损或外部灰尘等杂质不均匀黏附等而使转子的质量中心发生偏移, 引发了不平衡故障, 导致风机异常。

4 结论

本文分别采用时域分析、频域分析、小波分析对不平衡振动信号进行了研究。通过这三种方法的分析和对比, 有效地提取了不平衡故障的特征, 并对风机故障进行了诊断, 解决了故障问题, 以免给工矿企业带来不必要的损失。

参考文献

[1]张梅军.机械状态监测与故障诊断[M].北京:国防工业出版社, 2008.

快速傅里叶变化 篇6

1 FFT 算法的选择

对于N点有限长序列x( n) ,DFT( discrete fourier transform) 傅里叶变换对应的公式对为[4]:

式中,x( n) 代表n点的DFT序列,WnNk为DFT的旋转因子。FFT是DFT( 离散傅里叶变换) 的快速算法,它是根据DFT的奇、偶、虚、实特性演变而成。

传统的FFT算法包括卷积FFT算法和分裂FFT算法,卷积FFT算法是将DFT序列转换为卷积计算。分裂基算法是把有限长序列的N点DFT分解为每个长度为r的短序列进行DFT计算。分裂基算法相对卷积算法在逻辑控制上更加简单,这在硬件实现上占有很大优势。

FFT分裂基算法的表达式为N = rm,称为基r的N点FFT运算,其中m为FFT的级数。可见基二与基四、基八相比,FFT的序列长度N更灵活。

比较常用的基二和基四FFT算法,基四FFT的蝶形运算量为基二FFT的1 /4,级数为基二的一半,并且基四FFT的运算速度要比基二FFT块。但是,基四每个蝶形运算需要用到了三个复数乘法,而基二只需要一个复数乘法,从占用的FPGA资源方面考虑,基二算法更加节省资源。不仅如此,由于基四算法的复杂度较高,在一定程度上增加了器件连线长度和布线难度,使得芯片的连线延迟较高。综合这些因素,选择基二算法进行FFT硬件处理较为合理[5]。

基二FFT按时域抽 取法 ( DIT) 的蝶形运 算公式[6]:

式中,Am,Bm对应前一级的两个复数点,Am +1,Bm +1对应后一级两个复数点,WkN为旋转因子,该值可以预先计算得到。时域抽取法( DIT) 是在时域上分解N点DFT序列,输入为位倒序,输出为顺序。由上述公式可知,一次蝶形运算用到了1个复数乘法和2个复数加法。

2 1024点FFT硬件实现

硬件实现一般有顺序结构、并行结构、流水结构以及级联结构。其中顺序结构蝶形运算资源占用最少但速度最慢,并行结构最快但占用资源最大,流水结构虽然处理速度快但所占用资源仍然较大,而级联结构在资源占用与速度比例相对均衡。本文采用级联结构实现1024点FFT运算。

FFT的设计由时序控制单元,蝶形运算单元,地址管理单元,旋转因子存储单元及两组memory组成,如图1所示。

2. 1 旋转因子的产生与存储

旋转因子,既是X( k) 对x( n) 的离散傅里叶变换在频率处的采样值,也是Z平面单位圆上幅角在Z平面N等份后的第k点的值[7]。

因此,旋转因子 可以表示 为:

这里N为当前计算的点长度。第一级的旋转因子为W01024,W11024;第二级的旋转因子为W01024,W11024,W21024,W31024; 第九级的旋转因子为W01024,W11024,…,W11002234。这些旋转因子一一对应于Z变换域单位圆上,且各级旋转因子复用了前一级的旋转因子并使用了新的旋转因子,因此仅需存储第九级1 024个旋转因子即可。

将第九级产生的1,2,3,4象限的旋转因子的虚部实部分别存储到ROM中。而由上图可知,根据Z平面单位圆的对称性,第1、2、3、4象限的幅值是一样的,仅相位不同。第一象限的旋转因子分解到坐标轴后与其他象限分解后的值仅存在正负性差别,因此设计时,只需要存储第一象限的旋转因子即256个旋转因子即可,其他象限的旋转因子可通过对称性得到,进而节省所需的硬件资源。

2. 2 蝶形运算单元

蝶形运算是FFT运算的重要单元,该单元包括了乘法器电路,加法器,移位寄存器,饱和处理。如式( 3) 、式( 4) 所示,蝶形运算中需要将复数数据与旋转因子做复数乘法,及与其他数据加、减等运算。其中一次复数乘法如公式( 5) 所示,两组复数 ( A +Bj) 和 ( C + Dj) 相乘,包括4次实数乘法。而在式( 6) 、式( 7) 中对复数乘法做了优化[8]。

这种优化方式,使得原有的4次实数乘法转变为3次实数乘法,1次实数加法及2次实数减法。

硬件实现这种算法,只需提前知道C与D的值,得到C - D ,C + D即可。可以将上述公式的A,B对应蝶形运算输入数据的实部与虚部,C,D对应为旋转因子的实部与虚部,而由于FFT的旋转因子已知,即可通过查表的方式直接读出旋转因子的值并做加减运算,以式( 6) 、式( 7) 的方法完成蝶形运算复数乘法,节约乘法器资源。

本文对FFT数据采用定点化处理,定点化处理的优点在于实现简单,不需要额外的硬件资源。当单个蝶形运算单元的输入为N位时,输出最多为N + 2位,本文为了简化处理,将输入输出数据位宽都定为N + 2位。本文对每一级结果进行溢出判断,处理过程中对乘法结果进行了四舍五入,这样减小了截断误差。

2. 3 memory 存储单元

1 024点的FFT采用基二算法结构,需要10级FFT处理。本文采用8块256×28的双口RAM处理每级间数据,每4块RAM合成一组大小为1 024×28的memory,共合成两组,其中一组用于存储第0,2,4,6,8级的数据结果,另一组用于存储第1,3,5,7,9级的数据结果,以形成乒乓结构。如图2所示。

每一级的蝶形运算需要从对应的memory中的4块RAM里分别读取两个数据与旋转因子进行运算,每次并行处理八个数据,即4组蝶形运算,并在一个周期内完成,完成后将这8个点的结果写到另一组memory中,这样在128个周期后1 024个点都运算结束,如此往复,十级蝶形运算后,结束FFT。

2. 4 时序控制单元

时序控制单元包括状态机处理、数据控制逻辑。其中状态机主要处理每级蝶形运算间的状态跳转;数据控制逻辑主要负责memory数据读写控制以及不同memory间的乒乓操作控制。

由图3可知,第0级并行处理四个蝶形运算,每个蝶形运算单元处理同一块RAM里相邻地址的两个数据;

第1级的每个蝶形运算单元处理同一块RAM里间隔1个地址的两个数据;

第2级的每个蝶形运算单元处理同一块RAM里间隔3个地址的两个数据;

以此类推,第N( N < 8) 级时,则处理同一块RAM里间隔 ( 2N- 1) 个地址的两个数据;

当N≥8时,读取数据已超过了RAM的深度256,因此不能按照上述方法去处理。

当N = 8时,根据图3所示,从ram0和ram1的两个相同地址分别取两个数据,然后使用不同RAM但地址相同的两个数进行蝶形运算,ram2和ram3也按照此方式进行运算。

当N = 9时,从ram0和ram2的两个相同地址分别取两个数据,然后使用不同RAM但地址相同的两个数进行蝶形运算,ram1和ram3也按照此方式进行运算。

2. 5 地址生成单元

地址生成单元主要完成RAM的读写地址以及ROM的地址生成。地址生成单元包括计数器,截位拼接单元,块计数器,输出单元。计数器计算每一级所用周期数,每一级都为128 + 4个周期; 截位拼接单元用于通过计数器生成读取旋转因子的地址及读取memory的地址; 块计数器计算总共FFT的级数,范围为0 ~ 9。

需要注意的是,由于ROM中仅存有四分之一的旋转因子,因此需要根据需求对ROM中读出的旋转因子的值进行正负变换以得到符合要求的旋转因子。

3 时序仿真及结果

本文在ASIC下进行FFT的仿真及综合。将双口RAM及蝶形运算单元例化在FFT模块的测试代码中,双口RAM需要例化8块,进行ASIC仿真。以下两幅图为FFT仿真图。

先将输入数据写进memory2中,启动FFT的start信号,进行两组memory乒乓操作,图4( a) 即为FFT仿真结果。如图4( a) 所示,m从0级累计到9级,表示FFT内部控制memory数据读写的过程。

ASIC中的仿真结果显示,FFT每一级处理蝶形运算用到128个周期,则十级需要的时间是1 280周期,加上时序控制的状态转移30个周期,一共用1 310个周期。

ASIC的综合结果显示,整个FFT实现仅占用5万门资源。

本文用Xilinx公司的Virtex-7系列芯片做FFT的FPGA原型验证。Xilinx公司7系列的芯片工艺为28 nm,降低其成本和功耗,同时提高了FPGA的性能和容量。它通过整合6-LUT并与ARM公司合作规范AMBA总线,大幅提高芯片灵活性和片上系统性能。1024点的FFT数据,由于FPGA调试包括了DSP的控制,时钟频率仅开到50 MHz,仿真时间约为26. 2μs。如果在只有FFT硬件模块而没有DSP控制,时钟频率可以达到200 MHz,仿真时间仅有6. 55μs。

最后在matlab中做RTL硬件实现的误差分析。图5( b) 所示,real_error及image_error两个曲线为matlab结果与FPGA的仿真结果对应实部虚部的相对误差曲线,曲线在0附近波动,误差最大为5×10- 3,基本符合理论要求。图5( b) 中的matlab_real与dut_real曲线图分别对应matlab的FFT实部结果和RTL的FFT的实部结果,两条曲线基本一致,进一步说明RTL设计结果符合理论要求。

4 结论

通过对FFT一般算法的研究对比,提出了基2算法时域抽取法处理1 024点FFT的设计方案。通过了ASIC仿真、综合及FPGA的原型验证。该设计具有如下特点:

1) 四个蝶形运算并行处理,一个周期完成一次蝶形运算。

2) 用8块深度为256的双口memory读写FFT数据,并行8点处理FFT内部数据,共用1 320个周期。

3) 通过ASIC综合,硬件占用资源仅为5万门。

4) 用Xilinx的FPGA仿真测试,在50 MHz时钟频率下,需要时间为26. 4μs。

该FFT硬件实现是全RTL实现,已经通过前仿后仿及综合,可以作为IP核直接运用到实际工程中。

参考文献

[1] 解翀.基于FPGA的高速FFT处理器的研究.南京:南京航空航天大学,2010Xie C.Research on a high speed FFT processor with FPGA.Nanjing:Nanjing University of Aeronautics and Astronautics,2010

[2] 赵国亮.基于FPGA的1024点流水线结构FFT的算法研究与实现.西安:西安电子科技大学,2011Zhao G L.Research and implementation of 1024-point pipeline FFT algorithm based on FPGA.Xi’an:Xidian University,2011

[3] 李雪.基于FPGA的FFT处理器的设计与优化.哈尔滨:哈尔滨工业大学,2008Li X.The design and optimization of FFT processor on FPGA.Harbin:Harbin Institute of Technology,2008

[4] 谢彦林.可变点流水线结构FFT处理器的设计及FPGA实现.西安:西安电子科技大学,2007Xie Y L.The design of variable point pipeline FFT and implementation based on FPGA.Xi’an:Xidian University,2007

[5] 李小进,初建朋,赖宗声,等.高速基2FFT处理器的结构设计与FPGA实现.电路与系统学报,2005;10(5):49—53Li X J,Chu J P,Lai Z S,et al.FPGA implementation of a highspeed base-2 256-point FFT.Journal of Circuits and Systems(Natural Science),2005;10(5):49—53

[6] 王远模,赵宏钟,张军,等.用FPGA实现浮点FFT处理器的研究.国防科技大学学报,2004;26(6):61—64Wang Y M,Zhao H Z,Zhang J,et al.Realization of Floating-point FFT Processor with FPGA chip.Journal of National University of Defense Technology,2004;26(6):61—64

[7] 何方白,张德明,杨莉,等.数字信号处理.北京:高等教育出版社,2009:83—92He F B,Zhang D M,Yang L,et al.Digital signal processing.Beijing:Higher Education Press,2009:83—92