快速付里叶变换论文（共5篇）
快速付里叶变换论文 篇1
摘要：提出一种基-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个逆序数,即
为了便于比较,下面简单介绍生成法[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,即
后半部分中一个逆序数需要一次加法运算,即一次迭代共需要2i-1次加法,M此迭代共需要N-1次加法,即
所以,生成法一共需要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].西安工业学院学报.
[6]张学智,蔡晖.快速实现FFT的逆序方法[J].探测与控制学报,2001,23(2):62-64.
快速付里叶变换论文 篇2
关键词：模式识别,轮廓特征,2-范数,傅里叶,图形识别
1 引言
模式识别 (Pattern Recognition) :是指对表征事物或现象的各种形式的 (数值的、文字的和逻辑关系的) 信息进行处理和分析, 以对事物或现象进行描述、辨认、分类和解释的过程, 是信息科学和人工智能的重要组成部分。
近年来, 用于认知的人工智能 (模式识别) 得到了很大的发展和进步, 从环境的认知, 到人体的认知;从静态的认知到动态的认知。有关认知方面的人工智能, 例如:人脸检测、表情识别、车辆检测、语音识别等领域都已基本迈向成熟, 步入应用阶段。
图像识别属于模式识别的研究领域, 利用计算机视觉特征和机器学习的知识对目标对象进行识别分析。目前, 在各国研究学者的努力下得到了极大的发展, 图形识别过程中处理逻辑较简单, 数据来源基本来自图像或视频流, 以像素点为基本处理元素[1]。
其中, 轮廓特征在众多特征表示中可以说是最基本、最容易得到的, 获取轮廓特征的主要及首要步骤就是轮廓点的采集。在传统的处理方法中, 直接数据来源就是图片, 在二值化、膨胀、腐蚀等预处理后, 各轮廓点坐标其实就是轮廓像素点的横、纵坐标值。每一个像素点与其相邻像素点之间的距离都是相同的, 即等步长点, 且可采集的轮廓点数目也是有限的。但在软件的实际应用中, 直接数据来源却不适合以图片作为输入而是把图形轮廓拐点 (图形端点) 当作输入参数。两种直接原始数据各有优缺点。图片做源数据的方式, 能够在最大程度上保留图形的原始数据特征, 并且可以从中提取到更多特征信息, 如轮廓特征、区域特征、光流特征等等, 但是预处理步骤较多, 容易受到噪声点和冗余数据的影响。而直接以轮廓拐点做源数据的方式, 可以直接跳过图片预处理阶段, 减少图片中噪声点的影响。以此为代价的是, 我们要自己补全其它轮廓点坐标, 更重要的是, 轮廓坐标是实数轴上的, 相邻点之间不一定是等步长的, 因此, 我们要考虑的是:每两个点之间的取点步长要取多少?取多少轮廓点合适?……
2 处理流程
本文研究的主要工作是当以图形轮廓点作为直接数据来源时, 探讨更合适的特征处理方式。
识别的主要处理过程如图1所示。
3 主要工作
3.1 采点
在传统处理方式中, 采点时通常会按顺或逆时针的顺序把所有轮廓点全部提取, 然后通过要采集的点数确定采点步长 (相邻两个采集的特征点之间相隔的轮廓点个数) 依次取出要处理的特征点。
但像前面所提到的, 以轮廓点作为输入时, 首先要处理的就是考虑用来采点的轮廓点数要取多少 (补点数目) 。但是, 本实验中的实际点列处理是在实数轴上的, 每个轮廓点间可采点数目是无限的, 要还原所有轮廓点是不现实的 (严重影响了算法在实际应用中的效率) 。因此本文中的处理方式是将等步长的概念由间隔特征点个数相同转为间隔特征点间的长度相同。
(1) 、确定特征点数目n;
(2) 、计算图形周长l;
(3) 、步长steplength=l/n;
确定步长后的采点方式有两种:1、直接按步长、边线端点计算下一个特征点, 其中不会特意保留轮廓端点;2、计算方式和1相同, 不同的是, 在计算特征点坐标时, 如果轮廓端点P与起点的距离小等于steplength, 则保留P作为当前特征点。
算法保留轮廓端点的特征点采集
输入:图形轮廓端点
输出:特征点集
(1) 、Feature Points=[]
(2) 、Feature Points (1) =End Points (1)
(3) 、Start Point=End Points (1)
End Point=End Points (2)
(4) 、根据特征点数k计算步长steplength for I=2:k
(5) 、根据Start Point、End Point、steplength计算下一个特征点add Point
(6) 、将在Start Point和add Point之间的轮廓端点全部加入Feature Points
(7) 、判断add Point是否超过End Point, 是, 则将图形下一条边的端点作为新的Start Point和End Point;否, 则把add Point作为新的Start Point
End
(8) 、检查Feature Points的特征点数目是否超过k个, 若是, 则随机删除多余的非轮廓端点。
算法结束
下面几张图为上述两种补点方式采完点后的图形与原始图形之间的对比。
由上面几幅图可知, 第二种补点方式与第一种方式相比更能保存原始图形的局部特征, 能有效减少噪声点的影响。
3.2 特征表示
常用的特征表示有两种:静态全局形状特征和局部的动态特征。静态全局形状特征包括:边界描述 (链码、傅里叶描述子、离散近似、累加角函数、椭圆傅里叶描述子等) 和区域描述 (基本形状量度、矩、特性、重构) ;局部的运动特征主要有:光流、运动方向、轨迹、位置、速度。
本文要探讨两种特征表示:2-范数描述子、一维离散快速傅里叶描述子。
假设Feature Points中的特征点坐标表示如下:
求出质心坐标 (形心) :
求出质心坐标后, 再按 (3-3) 式分别计算各轮廓点到质心的距离, 计算公式采用欧氏距离:
那么该图形轮廓的边界-质心距可表示为:
3.2.12-范数描述子
2-范数对D做归一化处理以消除空间尺度和距离长度的影响, 2-范数公式
向量D即为对应图形处理后的2-范数特征描述子。
3.2.2一维离散快速傅里叶变换
傅里叶描述子 (Fourier Descriptor, 简称FD) , 是模式识别最常用的轮廓特征, 是物体形状边界曲线的傅立叶变换系数, 它是物体边界曲线信号的频域分析的结果[2]。
傅里叶描述子 (FD) 常用来表示封闭曲线的形状特征, 其基本思想是将目标轮廓曲线建模成一维序列, 对该序列进行一维的傅立叶变换, 从而获得一系列的傅立叶系数, 用这些系数对该目标轮廓进行描述。傅立叶描述子方法有一系列优点, 如:计算原理简单, 描述清晰, 具有由粗及精的特性等[3]。
由于傅立叶变换将序列的主要能量集中在了低频系数上, 因此, 傅立叶描述子的低频系数反映了轮廓曲线的整体形状。
一维函数f (x) (其中x=0, 1, 2, …M-1) 的傅里叶变换的离散形式为:
一维离散傅里叶变换虽然原理简单, 也比较容易实现, 但在实际应用中直接实现效率较低, 因此在实际工程中, 迫切需要一种能够快速计算离散傅里叶变换的高效算法, 快速傅里叶变换 (Fast Fourier Transform, 简称FFT) 由此而生, 其基本思路就是将较长的序列转换成相对短得多的序列来大大减少运算量[4]。
本文实验中直接使用MATLAB软件中自带的一维离散快速傅里叶变换函数fft进行计算。
由图4可知, 取低频30维部分作为特征描述即可。
4 实验结果分析
图5为采不同轮廓点时, 2-范数描述子和一维离散快速傅里叶描述子识别正确率的比较结果。由图5可知, 2-范数识别正确率在90%~92%之间浮动, 在2-范数处理时, 如果要取得更好的实验结果, 则轮廓点数不能太少 (通常在200左右做取舍) , 但随着轮廓点数的增加, 处理难度和算法效率都会受到影响。傅里叶描述子在95%上下浮动, 识别效果比2-范数更理想、更稳定。而且该特征描述子对参与计算的轮廓点数目要求不大, 从而也大大减少了计算消耗的时间。
图6为3.1节中两种取点方式实验结果对比, 其中, 由于第二种取点方式 (保留轮廓端点) 能更完整地保留原始图形的局部特征, 因此其识别效果比第一种取点方式 (不特意保留轮廓端点) 要稳定。若采用局部特征作为图形的表示特征, 那么这两种取点方式的优劣性会更明显一些, 因为局部特征点对全局特征表示的影响相对较小。
图7和图8分别是第一种采点方式和第二种采点方式做十折交叉验证平均结果的混淆矩阵。
十折交叉验证:将样本数据集随机分割为大小相等的十份, 每次取其中一份作为预测集, 剩下九份数据作为训练集, 计算当前识别正确率, 循环计算十次, 将十次预测正确率的平均值作为最后的结果。
由于La、Lb、Lc、Ld四种图形的特点限制, 其样本数目在总体样本集中的比例要远小于其它其余图形, 因此在做十折交叉验证时很容易识别成其它图形 (在实际应用中由于训练集相对实验时较完整, 因此不会出现这种情况) , 但第二种采点方式的识别结果仍然要好于第一种, 并且一维离散快速傅里叶描述子对其余图形的识别正确率比较理想 (训练集中样本数据要尽量发布均匀, 不能相差太多) 。
5 总结与展望
图形轮廓拐点 (图形端点) 作为直接输入参数时, 保留图形轮廓端点的一维离散快速傅里叶变换识别结果更加理想, 更为稳定, 能够很好地保留局部特征值, 同时也为以后增加新特征描述子的处理提供了更完整的数据源 (例如采用局部特征描述) 。
但是, 由于当前处理方法要自己模拟原始轮廓点列的排列情况, 因此只能处理轮廓为直线段或规则圆弧的对象, 在后续处理工作中可以利用曲线拟合方法模拟不规则图形轮廓表示图像, 再进行取点识别。
另外, 由于直接输入数据是图形轮廓端点, 除轮廓特征外, 很难得到其它特征表示, 尤其是局部特征的处理。可以考虑用自增的轮廓点模拟以图像流作为直接输入的情形, 从而获取部分局部特征, 用混合特征进行图形识别, 能取得更好的实验效果[5]。
在决策机制上, 可以利用简单集成学习方法对预测结果进行决策处理, 从而避免噪声数据, 增强识别的鲁棒性[6,7]。
参考文献
[1]黄静.基于形状特征的人体行为识别方法研究[D], 燕山大学, 2010.10~11
[2]王涛, 刘文印, 孙家广, 张宏江.傅里叶描述子识别物体的形状.计算机研究与发展.2002年.第39卷.1714~1719
[3]Dengsheng Zhang and Guojun Lu.Studyand evaluation of different Fourier methods for image retrieval.Image and Vision Computing.2005, 23:33-49.
[4]张铮, 王艳平, 薛桂香.数字图像处理与机器视觉.第5版.北京市崇文区夕照寺街14号:人民邮电出版社, 2012
[5]郭利, 姬晓飞, 李平, 曹江涛.基于混合特征的人体动作识别改进算法.[J].计算机应用研究, 2013)
[6]杨长盛, 陶亮, 曹振田, 汪世义.基于成对差异性度量的选择性集成方法[J].模式识别与人工智能, 2010, 第四期, 565~571
快速付里叶变换论文 篇3
关键词：心律失常,快速傅里叶变换,归一化相关系数,模板匹配
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.
快速付里叶变换论文 篇4
关键词：图像复原,快速全变差去卷积 (FTVd) ,快速傅里叶变换
1、引言
图像在形成、传输和记录的过程中, 由于很多原因使图像的质量变坏。典型的图像复原是根据图像退化的先验知识建立退化模型, 以此模型为基础, 采用各种逆退化处理方法进行复原。不失一般性, 给出图像复原的线性退化模型:
其中, u∈Rn2为原始图像, f∈Rn2为退化图像, K∈Rn2×n2为模糊退化算子, η∈Rn2为随机噪声。求解图像复原问题的方法很多, 常用的即求解一个最小二乘问题:为了求解的方便加入相应的正则项, 由Rudin等[1]人提出的全变差模型:
针对全变差模型 (2) 的解决[2], 仿真结果表明, FTVd算法的图像复原效果比较好。
2、全变差图像复原
对于复原模型 (2) , 应用优化理论中的变量分离与惩罚方法, 引入新的辅助量来代替梯度项Du, 产生约束的凸最小问题的新方程:
上式中的第二项称为惩罚项, β为惩罚因子。
符号说明:记D (1) , D (2) 分别为水平和垂直方向的一阶有限差分算子, D (1) u∈Rn2与D (2) u∈Rn2分别表示对图像u的水平与垂直方向的有限差分。记w= (w 1;w) ∈R2n2, 其中w1, w2∈Rn2分别是对D (1) u, D (2) u的一个近似, 让wi= ( (w 1) i; (w2) i) ∈R2是对Di u=[ (D (1) u) i, (D (2) u) i]∈R2的近似。
先固定u, (3) 式前两项关于iw可以分离, 因此 (3) 式关于w的最小值相当于求解
再固定w, 而u在 (3) 式是二次的, 因此u的最小值可以通过法方程给出:
在周期边界条件下, u, (D (1) ) T D (1) , (D (2) ) TD (2) 和KT K都是块循环矩阵, 因此 (5) 式的左边可以通过二维离散傅里叶变换来对角化。根据傅里叶变换的卷积理论[3], (5) 式可变为:
其中F表示快速傅里叶变换, *表示共轭转置, o表示分量方式乘积。每次迭代只有1w和2w是常量, 其余都是变量, 这些变量确定后, 就可以很快的计算出u。综上所述, 基于全变差的FTVd复原算法步骤如下:
(1) 输入f, tol>0, λ>0
1) 保存之前的迭代:up=u
2) 固定u, 根据 (4) 式计算出w
3) 固定w, 根据 (6) 式计算出u
(4) 输出:u
3、仿真结果
所有的实验都是在MATLABr2010a环境下实现的。图1给出了运动模糊的情况: (1) 是大小为的原始图像, (2) 是根据模型 (1) 获得的256×256退化图像, SNR为7.2d B。模糊算子K是运动核函数, 随机噪声为高斯白噪声, 方差为10-7; (3) 是基于全变差的FTVd算法, 设定β的初值为1, 终值βmax=220, β的增长率为2, 内部循环迭代容许的误差tol=5.0×10-4;经过 (3) 复原后得到的图像SN R为22.19d B, CPU所用时间为0.5秒, 图像的质量显然得到了很大的提高。
4、总结
基于全变差的图像复原是现今图像复原的热点之一, 对于本文中的FTVd算法还有待于改进, 算法中要求不断地增加惩罚因子β的值, 随之问题的条件数就会变大, 会使得问题越难求解, 同时会带来数值不稳定问题。这些问题我们将在以后的工作解决。
参考文献
[1]L.Rudin and S.Osher.Nonlinear total variation based noise removal algorithms, PhysD, (1992) , pp.259-268.
[2]L.Rudin and S.Osher.Total variation based image restoration with free local constraints.Proc.1st IEEE ICIP, 1:31-35, 1994.
[3]C.Van Loan.Computational Frameworks for the Fast Fourier Transform, volume 10 of Frontier in Applied Mathematics.SIAM, 1992.
快速付里叶变换论文 篇5
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]。
因此,旋转因子 可以表示 为:
将第九级产生的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