CN100576200C - 一种基于多抽样的分数阶傅立叶变换实现方法 - Google Patents
一种基于多抽样的分数阶傅立叶变换实现方法 Download PDFInfo
- Publication number
- CN100576200C CN100576200C CN200810106727A CN200810106727A CN100576200C CN 100576200 C CN100576200 C CN 100576200C CN 200810106727 A CN200810106727 A CN 200810106727A CN 200810106727 A CN200810106727 A CN 200810106727A CN 100576200 C CN100576200 C CN 100576200C
- Authority
- CN
- China
- Prior art keywords
- delta
- fpga
- chirpb
- chirpa
- signal processing
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Complex Calculations (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明属于信号处理领域,具体涉及一种基于FPGA的分数阶傅立叶变换信号处理实现方法,提高在硬件平台实时实现的计算效率。其基本原理是在土耳其人Ozaktas离散采样型算法基础上,结合FPGA平台并行结构和流水处理特点,采用多抽样信号处理理论,将原始算法分解成两路并行的多相等效结构,去除了冗余运算,优化了算法流程,解决了原有方法实现时计算速度慢,时延大,不利于实时处理等问题,相对原始算法流程计算量更小,效率更高,并且由于是并行结构,更适合FPGA等硬件实现;另外,输入序列预处理的插值滤波部分选用了加窗半带滤波器,不仅解决了低阶次时出现的边缘振荡问题,而且由于滤波器长度较短,有利于FPGA平台流水处理,提高实现时的计算效率。
Description
所属技术领域
本发明属于信号处理领域,具体涉及一种基于FPGA(Field ProgrammableGate Array,现场可编程门阵列)的分数阶傅立叶变换信号处理高效实现方法,提高在硬件平台实时实现的计算效率。
背景技术
分数阶傅立叶变换最初在光学领域具有广泛应用,1993年Almeida把分数阶傅立叶变换解释为信号在时频平面的旋转,是经典傅立叶变换的推广;1996年Ozaktas提出了一种与FFT计算速度相当的离散采样型算法后,分数阶傅立叶变换才开始在信号处理领域得到应用。分数阶傅立叶变换可以看成是一种统一的时频变换,同时反映了信号在时、频域的信息,与常用二次型时频分布不同的是它用单一变量来表示时频信息,且没有交叉项困扰,与传统傅立叶变换(其实是分数阶傅立叶变换的一个特例)相比,它适于处理非平稳信号,尤其是chirp类信号,且多了一个自由参量(变换阶数a),因此分数阶傅立叶变换在某些条件下往往能够得到传统时频分布或傅立叶变换所得不到的效果,而且由于它具有比较成熟的快速离散算法,因此在得到更好效果的同时并不需要付出太多的计算代价。目前信号处理领域对分数阶傅立叶变换的应用有如下4种方式:
(1)分数阶傅立叶变换是一种统一的时频变换,随着阶数从0连续增长到1,分数阶傅立叶变换展示出信号从时域逐步变化到频域的所有变化特征,可以为信号的时频分析提供更大的选择余地;最直接的利用方式就是将传统时、频域的应用推广到分数阶傅立叶域以获得某些性能上的改善,如:分数阶傅立叶域的滤波等。
(2)分数阶傅立叶变换可以理解为chirp基分解,因此,它十分适合处理chirp类信号,而chirp类信号在雷达、通信、声纳以及自然界中经常遇到。分数阶傅立叶变换可以处理这些领域中的波束形成、目标识别以及扫频干扰抑制等等问题。
(3)分数阶傅立叶变换是对时频平面的旋转,利用这一点可以建立起分数阶傅立叶变换与时频分析工具的关系,既可以用来估计瞬时频率、恢复相位信息,又可以用来设计新的时频分析工具,如:TTFT、以高斯函数的分数阶傅立叶变换为基函数的信号扩展方法等。
(4)与傅立叶变换相比,分数阶傅立叶变换多一个自由参数,因此在某些应用场合能够得到更好的效果,如:数字水印和图像加密。
分数阶傅里叶变换定义式如下:
Ba(u,x)≡Aφ exp[iπ(u2cotφ-2ux cscφ+x2 cotφ)],
其中Aφ≡{exp(-iπ sgn(sinφ)/4+iφ/2)}/|sinφ|,φ≡aπ/2。
重新整理上式可以得到:
分数阶傅立叶变换的离散算法有很多种,其中土耳其人Ozaktas提出的快速算法是目前主流算法的一种,其基本原理是将经过量纲归一化以后的信号f(x)的分数阶傅立叶变换分解成以下3个简单步骤,并对每个步骤进行离散化;
1.首先用chirp信号调制信号f(x)得到g(x):
2.调制信号g(x)与另一个chirp信号卷积:
3.用chirp信号调制卷积后的信号:
关于量纲归一化可以参考文献“分数阶傅立叶变换数值计算中的量纲归一化研究”,发表于《北京理工大学学报》2005年第四期,作者是赵兴浩、邓兵、陶然)。
Ozaktas提出的快速算法的具体实现过程为:
其中N=Δ2。将(2)式代入(1)式化简得到:
然后以1/2Δ为采样间隔,在[-Δ/2,Δ/2]范围内对分数阶域进行离散化,令 则:
其中-N≤m≤N,进一步整理上式得到:
假设原始信号序列长度为N,归一化采样率 采样间隔为1/Δ,而以上3个步骤中g(x),Ga(u)的采样间隔都是1/2Δ,因此在进行第一步chirp调制之前需要对原始信号样本序列进行2倍内插,具体做法是先将信号序列零值内插,然后通过一个低通滤波器h(n)平滑,经过零值插值的结果与2N长度chirp序列相乘,再与4N长度序列作卷积,卷积结果截取2N点再与第一个chirp序列及一个常系数相乘,最后还要进行2倍抽取以便得到采样间隔为1/Δ的变换结果序列类似FFT,总的变换过程是N输入,N输出的。
按照上述算法流程的经典实现方法存在以下两个缺点:
(1)运算过程速度慢,时延大。一个N输入,N输出的变换过程中间需要经过一次4N长度的卷积,也就是3次4N长度的FFT或IFFT变换,对硬件平台资源要求较高,不利于FPGA等实时实现。
(2)目前应用较多的零值内插方法主要有以下两种,各自有不同的缺点:一种是sinc函数序列作低通滤波器平滑零值内插序列,一般这个滤波器长度为输入序列长度的4倍,效果好,但是计算量大,不利于FPGA等硬件实现,而且阶次较低时,变换结果的边缘出现振荡;另一种是用拉格朗日插值方法,这种方法计算速度快,但效果不好,尤其是在计算结果的边缘误差大。A.Bultheel对这两种方法进行了比较,具体请看文献“Computation of the Fractional Fouriertransform”,发表于《Applied and Computational Harmonic Analysis》2004年第16期。
发明内容
为了解决上述经典实现方法中计算速度慢,时延大,不利于实时处理等问题,本发明提出了一种基于FPGA的分数阶傅立叶变换信号处理高效实现方法,相对原始算法流程计算量更小,效率更高,并且由于是并行结构,更适合FPGA等硬件实现;另外,本发明在输入序列预处理的插值滤波部分选用了更有效的加窗半带滤波器,不仅解决了低阶次时出现的边缘振荡问题,而且由于滤波器长度较短,有利于FPGA平台流水处理,提高实现时的计算效率。
本发明方法的基本原理是在Ozaktas提出的快速算法基础上,结合FPGA平台并行结构和流水处理特点,采用多抽样信号处理理论,将原始算法分解成两路并行的多相等效结构,去除了冗余运算,优化了算法流程,提高了计算效率。
具体发明内容如下:
一种基于FPGA的分数阶傅立叶变换信号处理实现方法,其特征在于,包含以下5个步骤,每个步骤包括两个并行的信号处理支路:
①在FPGA中,将输入信号序列f(n/Δ)通过插值滤波器h(n)的多相分量ri(k)进行滤波,得到序列fi(m/Δ),其中i=0,1表示支路数。|n|≤(N-1)/2, 表示输入信号序列归一化的采样间隔,N表示输入信号序列长度,为奇数,其中:
h(n)=2·ω(n)·sin(πn/2)/(πn),|n|≤L,2L+1为滤波器长度;
ω(n)=0.5·[1+cos(nπ/L)],|n|≤L;
r0(n)=h(2n+1),r1(n)=h(2n)=δ(n);
②用FPGA中的Chirp DDS产生chirp信号chirpA(n/(2Δ))和chirpB(n/(2Δ)),用chip信号chirpA(n/(2Δ))的两个多相分量li(k/Δ)分别与步骤①相应支路的滤波结果fi(m/Δ)相乘得到gi(k/Δ),其中
③gi(k/Δ)与另一chirp信号chirpB(n/(2Δ))的两个多相分量ei(k/Δ)分别作卷积,卷积采用FPGA中的FFT模块实现,其中
④将步骤③中的两路卷积结果分别截取N点;
⑤将步骤④得到的两个支路的N点卷积结果相加,相加结果再与步骤②中的序列l0(k/Δ)和系数Aφ/(2Δ)相乘,此处相乘的步骤采用FPGA中的IFFT模块实现,输出为最后结果其中,|m|≤(N-1)/2,Aφ≡exp(-iπsgn(sinφ)/4+iφ/2)/|sinφ|。
下面对上述发明内容的推导过程进行简要说明:
假设输入信号序列为:
归一化采样率为 经过插值以后,采样率变为2Δ,新序列长度变为2N-1:
设抗混叠插值滤波器为h(n)。
则零值内插后序列为:
抗混叠滤波后得到:
其中-(N-1)≤n≤N-1,-(N-1)/2≤k≤(N-1)/2。
设
r0(n)=h(2n+1)
r1(n)=h(2n)
当n=2m时:
当n=2m+1时:
又
其中φ≡aπ/2,a为变换阶次。上面两个chrip序列的多相分量分别如下:
则令
根据(5)式有
其中|m|≤(N-1)/2,Aφ≡exp(-iπsgn(sinφ)/4+iφ/2)/|sinφ|。(7)式即可按前面介绍的5个步骤实现整个计算过程,整个过程具有两个支路并行运算的特点。
需要说明的是,在上述推导过程中并没有限定插值滤波器的类型,即本推导过程对任意一种插值滤波器都适用。
本发明提出的一种基于FPGA的分数阶傅立叶变换信号处理实现方法,其有益效果在于:
(1)本发明提出的实现方法相比于原始算法经典实现方法运算量更小。设插值滤波器长度为4N,那么两种实现方法的运算量如表1所示,从表1中可以看出多抽样高效实现方法每一支路的运算量均小于原始算法运算量的一半;
(2)本发明提出的实现方法在结构上具有两路并行的特点,实现时可以采用并行结构来提高计算速度,节省计算时间,尤其适合于FPGA等硬件实现。具体可参见“具体实施方式”中的FPGA实现方法实例;
(3)本发明提出的实现方法选用低阶次的加窗半带滤波器作为插值滤波器,节省了插值滤波环节一半的计算量,并且加窗可以有效防止在变换阶次较低时出现的振荡现象,具体见附图2中的比较。从图中可以看出插值滤波器在阶次较低,不加窗情况下,分数阶变换结果出现了振荡,误差大,而加窗后,则消除了振荡。
表1本发明提出的实现方法和经典实现方法运算量的比较
运算量 | 原始算法经典实现 | 多抽样实现单一支路 |
复乘数 | 12N log<sub>2</sub> N+36N | (11N log<sub>2</sub> N+22N)/2 |
复加数 | 24N log<sub>2</sub> N+48N | (22N log<sub>2</sub> N+23N)/2 |
附图说明
图1为多抽样高效实现方法流程图;
图2为Hanning窗半带滤波器与其它插值滤波器对分数阶傅立叶变换结果影响的比较;
图3是一个输入信号长度为1023点的本发明方法的FPGA实现结构框图;
具体实施方式
下面结合附图及FPGA实施例对发明内容做详细说明。
本发明涉及到一种基于FPGA的分数阶傅立叶变换信号处理实现方法,其原理见式(7),实现方法的算法流程图如图1所示,整个流程分解成以下5个步骤完成运算:
⑥在FPGA中,将输入信号序列f(n/Δ)通过插值滤波器h(n)的多相分量ri(k)进行滤波,得到序列fi(m/Δ),其中i=0,1表示支路数。|n|≤(N-1)/2, 表示输入信号序列归一化的采样间隔,N表示输入信号序列长度,为奇数,其中
h(n)=2·ω(n)·sin(πn/2)/(πn),|n|≤L,2L+1为滤波器长度;
ω(n)=0.5·[1+cos(nπ/L)],|n|≤L;
r0(n)=h(2n+1),r1(n)=h(2n)=δ(n);
⑦用FPGA中的Chirp DDS产生chirp信号chirpA(n/(2Δ))和chirpB(n/(2Δ)),用chirp信号chirpA(n/(2Δ))的两个多相分量li(k/Δ)分别与步骤①相应支路的滤波结果fi(m/Δ)相乘得到gi(k/Δ),其中
⑧gi(k/Δ)与另一chirp信号chirpB(n/(2Δ))的两个多相分量ei(k/Δ)分别作卷积,卷积采用FPGA中的FFT模块实现,其中
⑨将步骤③中的两路卷积结果分别截取N点;
⑩将步骤④得到的两个支路的N点卷积结果相加,相加结果再与步骤②中的序列l0(k/Δ)和系数Aφ/(2Δ)相乘,此处相乘的步骤采用FPGA中的IFFT模块实现,输出为最后结果其中,|m|≤(N-1)/2,Aφ≡exp(-iπsgn(sinφ)/4+iφ/2)/|sinφ|。
下面结合上述5个步骤给出一个该算法用于FPGA实现的例子,输入信号的长度设为1023点,插值滤波器是123阶的Hanning窗半带FIR滤波器。图3是本例的FPGA实现原理框图,图中省略了控制模块,由以下几个基本模块组成:插值滤波模块;ChirpDDS模块;FFT和IFFT模块;输出单元;控制模块。
(1)插值滤波模块是一个123阶的Hanning窗半带FIR滤波器h(n),由于采用多抽样结构,r1(n)=δ(n),因此第二支路的插值滤波模块可以省略。再考虑系数的对称性,第一支路只需一个31阶的FIR滤波器就可实现,其中
(2)ChirpDDS模块采用传统的Chirp DDS设计方法,由相位累加单元和频率累加单元再加上控制单元和正弦查找表构成。主要功能是输出公式(6)中对应的几个chirp序列:chirpA(n/2Δ),chirpB(n/2Δ)。序列chirpA的两个多相分量l0(k/Δ),l1(k/Δ)用于步骤②两个支路的调制。序列chirpB的两个多相分量e0(k/Δ),e1(k/Δ)用于步骤③的卷积运算,其中
(3)FFT模块和IFFT模块完成步骤③的卷积运算,两个支路各一个FFT模块,最后合用一个IFFT模块。3个模块进行FFT变换的长度是由控制器控制的,最大长度都是输入序列长度的2倍,即2N≈2048;
(4)输出单元由一个复乘法器和一个缓存序列l0(k/Δ)及常数的ram组成。对应步骤④和⑤。从ram中读取常序列与IFFT输出也就是卷积输出相乘得到最终的1023点分数阶傅立叶变换结果;
(5)控制单元是一个复杂的状态机,它的主要功能如下:根据变换阶次的大小控制ChirpDDS输出序列的调频率和起始频率;控制FFT和IFFT模块的变换长度;控制ram的读写使能信号、读写地址以及其他控制信号。
假设现在需要对一帧N=1023点数据的进行连续100阶次分数阶傅立叶变换,整个模块在Xilinx xc2vp20 FPGA芯片中实现,工作在100M时钟,那么整个计算过程消耗的时间大约为4.1ms,这要比传统方法的FPGA实现的9.5ms快一倍左右。
以上所述的具体描述,对发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限定本发明的保护范围,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (1)
1.一种基于FPGA的分数阶傅立叶变换信号处理实现方法,其特征在于,包含以下5个步骤,每个步骤包括两个并行的信号处理支路:
①在FPGA中,将输入信号序列f(n/Δ)通过插值滤波器h(n)的多相分量ri(k)进行滤波,得到序列fi(m/Δ),其中i=0,1表示支路数。|n|≤(N-1)/2, 表示输入信号序列归一化的采样间隔,N表示输入信号序列长度,为奇数,其中;
h(n)=2·w(n)·sin(πn/2)/(πn),|n|≤L,2L+1为滤波器长度;
w(n)=0.5·[1+cos(nπ/L)],|n|≤L;
r0(n)=h(2n+1),r1(n)=h(2n)=δ(n);
②用FPGA中的Chirp DDS产生chirp信号chirpA(n/(2Δ))和chirpB(n/(2Δ)),用chirp信号chirpA(n/(2Δ))的两个多相分量li(k/Δ)分别与步骤①相应支路的滤波结果fi(m/Δ)相乘得到gi(k/Δ),其中
③gi(k/Δ)与另一chirp信号chirpB(n/(2Δ))的两个多相分量ei(k/Δ)分别作卷积,卷积采用FPGA中的FFT模块实现,其中
④将步骤③中的两路卷积结果分别截取N点;
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN200810106727A CN100576200C (zh) | 2008-05-15 | 2008-05-15 | 一种基于多抽样的分数阶傅立叶变换实现方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN200810106727A CN100576200C (zh) | 2008-05-15 | 2008-05-15 | 一种基于多抽样的分数阶傅立叶变换实现方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN101303689A CN101303689A (zh) | 2008-11-12 |
CN100576200C true CN100576200C (zh) | 2009-12-30 |
Family
ID=40113597
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN200810106727A Expired - Fee Related CN100576200C (zh) | 2008-05-15 | 2008-05-15 | 一种基于多抽样的分数阶傅立叶变换实现方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN100576200C (zh) |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101645865B (zh) * | 2009-09-15 | 2011-12-28 | 哈尔滨工业大学 | 基于分数阶傅立叶变换的信道中多径时延和多普勒频移估计方法及实现该方法的系统 |
CN103067329B (zh) * | 2012-12-25 | 2015-11-11 | 桂林电子科技大学 | 一种chirp实信号的快速傅里叶变换多带合成和分离方法 |
CN103685128B (zh) * | 2013-12-27 | 2017-04-12 | 湖北民族学院科技学院 | 应用于ofdm发射机的ifft处理器和ifft实现方法 |
IE87469B1 (en) | 2016-10-06 | 2024-01-03 | Google Llc | Image processing neural networks with separable convolutional layers |
CN107168927B (zh) * | 2017-04-26 | 2020-04-21 | 北京理工大学 | 一种基于流水反馈滤波结构的稀疏傅里叶变换实现方法 |
CN109683478A (zh) * | 2018-12-21 | 2019-04-26 | 南京埃斯顿机器人工程有限公司 | 柔性关节机械臂分数阶滑模优化控制方法 |
-
2008
- 2008-05-15 CN CN200810106727A patent/CN100576200C/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN101303689A (zh) | 2008-11-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN100576200C (zh) | 一种基于多抽样的分数阶傅立叶变换实现方法 | |
CN103646011B (zh) | 一种基于线性调频z变换的信号频谱细化方法 | |
CN102999473A (zh) | 一种线性调频信号的检测与参数估计方法 | |
Ali Khan et al. | Sparsity-aware adaptive directional time–frequency distribution for source localization | |
Tao et al. | An efficient FPGA-based implementation of fractional Fourier transform algorithm | |
CN106842144A (zh) | 并行多相结构数字脉压方法 | |
CN104052494A (zh) | 面向频域稀疏信号的信号重构方法 | |
Borisagar et al. | Fourier transform, short-time fourier transform, and wavelet transform | |
Liu et al. | Optimal target function for the fractional Fourier transform of LFM signals | |
Alsteris et al. | Iterative reconstruction of speech from short-time Fourier transform phase and magnitude spectra | |
CN107632963A (zh) | 长度为复合数的基于fft的电力系统采样信号希尔伯特变换方法 | |
Crowley et al. | How hard is the euro area core? An evaluation of growth cycles using wavelet analysis | |
Wei et al. | Time–frequency analysis method based on affine Fourier transform and Gabor transform | |
CN101741801B (zh) | 一种32路并行数据dft的实现结构 | |
Xiong et al. | Fractional domain singularity power spectrum | |
Morabito et al. | Improved computational performance for distributed passive radar processing through channelised data | |
Abel et al. | Spectral analysis of atmospheric radar echoes using a non-stationary approach | |
Yu-long et al. | Modification and digital implementation of fam algorithm for spectral correlation | |
CN105656451A (zh) | 一种基于频域处理的扩频信号匹配滤波系统及方法 | |
Raju et al. | Fast implementation of sparse iterative covariance-based estimation for processing MST radar data | |
Zhen et al. | Spatial Spectrum Estimation for Wideband Signals by Sparse Reconstruction in Continuous Domain | |
CN104820581B (zh) | 一种fft和ifft逆序数表的并行处理方法 | |
Zhang et al. | Algorithm for finding the best linear combination window under the sparse solution constraint of discrete gabor transform coefficients | |
CN103810146B (zh) | 一种逆序输入顺序输出的fft结构设计方法 | |
CN101977032A (zh) | 应用于全数字b型超声诊断仪中的动态滤波器 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
C17 | Cessation of patent right | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20091230 Termination date: 20130515 |