首页> 中国专利> 一种基于多抽样的分数阶傅立叶变换实现方法

一种基于多抽样的分数阶傅立叶变换实现方法

摘要

本发明涉及一种基于多抽样的分数阶傅立叶变换实现方法,属于信号处理领域,用于提高实时实现的计算效率。该方法在离散采样型算法基础上,采用多抽样信号处理理论,将原始算法分解成两路并行的多相等效结构,去除了冗余运算,优化了算法流程,解决了Ozaktas离散采样型算法实现时计算速度慢,时延大,不利于实时处理等问题,相对原始算法流程计算量更小,效率更高,并且由于是并行结构,更适合硬件实现;另外,本发明在输入序列预处理的插值滤波部分选用了更有效的加窗半带滤波器,不仅解决了低阶次时出现的边缘振荡问题,而且由于滤波器长度较短,有利于流水处理,提高实现时的计算效率。

著录项

  • 公开/公告号CN101303689A

    专利类型发明专利

  • 公开/公告日2008-11-12

    原文格式PDF

  • 申请/专利权人 北京理工大学;

    申请/专利号CN200810106727.2

  • 发明设计人 陶然;梁广平;黄克武;单涛;

    申请日2008-05-15

  • 分类号G06F17/14(20060101);

  • 代理机构11120 北京理工大学专利中心;

  • 代理人张利萍

  • 地址 100081 北京市海淀区中关村南大街5号

  • 入库时间 2023-12-17 21:06:40

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2014-07-09

    未缴年费专利权终止 IPC(主分类):G06F17/14 授权公告日:20091230 终止日期:20130515 申请日:20080515

    专利权的终止

  • 2009-12-30

    授权

    授权

  • 2009-01-07

    实质审查的生效

    实质审查的生效

  • 2008-11-12

    公开

    公开

说明书

技术领域

本发明属于信号处理领域,具体涉及一种基于多抽样的分数阶傅立叶变换高效实现方法,提高实时实现的计算效率。

背景技术

分数阶傅立叶变换最初在光学领域具有广泛应用,1993年Almeida把分数阶傅立叶变换解释为信号在时频平面的旋转,是经典傅立叶变换的推广;1996年Ozaktas提出了一种与FFT计算速度相当的离散采样型算法后,分数阶傅立叶变换才开始在信号处理领域得到应用。分数阶傅立叶变换可以看成是一种统一的时频变换,同时反映了信号在时、频域的信息,与常用二次型时频分布不同的是它用单一变量来表示时频信息,且没有交叉项困扰,与传统傅立叶变换(其实是分数阶傅立叶变换的一个特例)相比,它适于处理非平稳信号,尤其是chirp类信号,且多了一个自由参量(变换阶数a),因此分数阶傅立叶变换在某些条件下往往能够得到传统时频分布或傅立叶变换所得不到的效果,而且由于它具有比较成熟的快速离散算法,因此在得到更好效果的同时并不需要付出太多的计算代价。目前信号处理领域对分数阶傅立叶变换的应用有如下4种方式:

(1)分数阶傅立叶变换是一种统一的时频变换,随着阶数从0连续增长到1,分数阶傅立叶变换展示出信号从时域逐步变化到频域的所有变化特征,可以为信号的时频分析提供更大的选择余地;最直接的利用方式就是将传统时、频域的应用推广到分数阶傅立叶域以获得某些性能上的改善,如:分数阶傅立叶域的滤波等。

(2)分数阶傅立叶变换可以理解为chirp基分解,因此,它十分适合处理chirp类信号,而chirp类信号在雷达、通信、声纳以及自然界中经常遇到。分数阶傅立叶变换可以处理这些领域中的波束形成、目标识别以及扫频干扰抑制等等问题。

(3)分数阶傅立叶变换是对时频平面的旋转,利用这一点可以建立起分数阶傅立叶变换与时频分析工具的关系,既可以用来估计瞬时频率、恢复相位信息,又可以用来设计新的时频分析工具,如:TTFT、以高斯函数的分数阶傅立叶变换为基函数的信号扩展方法等。

(4)与傅立叶变换相比,分数阶傅立叶变换多一个自由参数,因此在某些应用场合能够得到更好的效果,如:数字水印和图像加密。

分数阶傅里叶变换定义式如下:

Fa[f(u)]{Faf}(u)-Ba(u,x)f(x)dx,

Ba(u,x)≡Aφexp[iπ(u2cotφ-2ux cscφ+x2 cotφ)],

其中Aφ≡{exp(-iπsgn(sinφ)/4+iφ/2)}/|sinφ|,φ≡aπ/2。

重新整理上式可以得到:

{Faf}(u)=Aφeu2cotφ-e-j2πuxcscφ[ex2cotφf(x)]dx

分数阶傅立叶变换的离散算法有很多种,其中土耳其人Ozaktas提出的快速算法是目前主流算法的一种,其基本原理是将经过量纲归一化以后的信号f(x)的分数阶傅立叶变换分解成以下3个简单步骤,并对每个步骤进行离散化:

1.首先用chirp信号调制信号f(x)得到g(x):

g(x)=f(x)exp(-x2tanφ2).

2.调制信号g(x)与另一个chirp信号卷积:

Ga(u)=Aφ-exp[(u-x)2cscφ]g(x)dx

3.用chirp信号调制卷积后的信号:

{Faf}(u)=Ga(u)exp(-u2tanφ2)

关于量纲归一化可以参考文献“分数阶傅立叶变换数值计算中的量纲归一化研究”,发表于《北京理工大学学报》2005年第四期,作者是赵兴浩、邓兵、陶然)。

Ozaktas提出的快速算法的具体实现过程为:

假定信号f(x)的Wigner-Ville分布限定在原点为中心,直径为Δ的圆内。则当分数阶变换阶次a满足0.5≤|a|≤1.5时,的最大带宽限定在Δ内,以1/2Δ为采样间隔,以香农公式重建可得:

ex2cotφf(x)=Σn=-NNe(n2Δ)2cotφf(n2Δ)sinc(2Δ(x-n2Δ))---(2)

其中N=Δ2。将(2)式代入(1)式化简得到:

{Faf}(u)=Aφ2Δeu2cotφΣn=-NNe-j2πun2Δcscφe(n2Δ)2cotφf(n2Δ)---(3)

然后以1/2Δ为采样间隔,在[-Δ/2,Δ/2]范围内对分数阶域进行离散化,令u=m2Δ,则:

{Faf}(m2Δ)=Aφ2Δe(m2Δ)2cotφΣn=-NNe-j2πmn(2Δ)2cscφe(n2Δ)2cotφf(n2Δ)---(4)

其中-N≤m≤N,进一步整理上式得到:

{Faf}(m2Δ)=Aφ2Δe-(m2Δ)2tanφ2Σn=-NNejπ(m-n2Δ)2cscφe-(n2Δ)2tanφ2f(n2Δ)]---(5)

假设原始信号序列长度为N,归一化采样率Δ=N,采样间隔为1/Δ,而以上3个步骤中g(x),Ga(u)的采样间隔都是1/2Δ,因此在进行第一步chirp调制之前需要对原始信号样本序列进行2倍内插,具体做法是先将信号序列零值内插,然后通过一个低通滤波器h(n)平滑,经过零值插值的结果与2N长度chirp序列相乘,再与4N长度序列作卷积,卷积结果截取2N点再与第一个chirp序列及一个常系数相乘,最后还要进行2倍抽取以便得到采样间隔为1/Δ的变换结果序列{Faf(m/Δ),类似FFT,总的变换过程是N输入,N输出的。

按照上述算法流程的经典实现方法存在以下两个缺点:

(1)运算过程速度慢,时延大。一个N输入,N输出的变换过程中间需要经过一次4N长度的卷积,也就是3次4N长度的FFT或IFFT变换。

(2)目前应用较多的零值内插方法主要有以下两种,各自有不同的缺点:一种是sinc函数序列作低通滤波器平滑零值内插序列,一般这个滤波器长度为输入序列长度的4倍,效果好,但是计算量大,不利于硬件实现,而且阶次较低时,变换结果的边缘出现振荡;另一种是用拉格朗日插值方法,这种方法计算速度快,但效果不好,尤其是在计算结果的边缘误差大。A.Bultheel对这两种方法进行了比较,具体请看文献“Computation of the Fractional Fourier transform”,发表于《Applied and Computational Harmonic Analysis》2004年第16期。

发明内容

本发明的目的是为了解决上述经典实现方法中计算速度慢,时延大,不利于实时处理等问题,提出了一种基于多抽样的分数阶傅立叶变换实现方法,该方法相对原始算法流程计算量更小,效率更高,并且由于是并行结构,更适合硬件实现;另外,本发明在输入序列预处理的插值滤波部分选用了更有效的加窗半带滤波器,不仅解决了低阶次时出现的边缘振荡问题,而且由于滤波器长度较短,有利于流水处理,提高实现时的计算效率。

本发明方法的基本原理是在Ozaktas提出的快速算法基础上,采用多抽样信号处理理论,将原始算法分解成两路并行的多相等效结构,去除了冗余运算,优化了算法流程,提高了计算效率。

本发明的目的是通过下述技术方案实现的。

本发明的一种基于多抽样的分数阶傅立叶变换实现方法,包含以下5个步骤,每个步骤包括两个并行运算的支路:

①将输入信号序列f(n/Δ)通过插值滤波器h(n)的多相分量ri(k)进行滤波,得到序列fi(m/Δ),其中i=0,1表示支路数。|n|≤(N-1)/2,Δ=N表示输入信号序列归一化的采样间隔,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);

②用chirp信号chirpA(n/(2Δ))的两个多相分量li(k/Δ)分别与步骤①相应支路的滤波结果fi(m/Δ)相乘得到gi(k/Δ),其中:

chirpA(n2Δ)=exp(-jπ4(Δ)2n2tanφ2),|n|≤N-1,φ≡aπ/2,a为变换阶次;

l0(nΔ)=chirpA(2n2Δ),l1(nΔ)=chirpA(2n+12Δ);

③gi(k/Δ)与另一chirp信号chirpB(n/(2Δ))的两个多相分量ei(k/Δ)分别作卷积,其中:

chirpB(n2Δ)=exp(jπ4(Δ)2n2cscφ),|n|≤2N-1;

e0(nΔ)=chirpB(2n2Δ),e1(nΔ)=chirpB(2n+12Δ);

④将步骤③中的两路卷积结果分别截取N点;

⑤将步骤④得到的两个支路的N点卷积结果相加,相加结果再与步骤②中的序列l0(k/Δ)和系数Aφ/(2Δ)相乘,输出最后结果{Faf}(m/Δ),其中|m|≤(N-1)/2,Aφ≡exp(-iπsgn(sinφ)/4+iφ/2)/|sinφ|。

下面对上述发明内容的推导过程进行简要说明:

假设输入信号序列为:

|n|N-12,N表示序列长度为奇数

归一化采样率为Δ=N,经过插值以后,采样率变为2Δ,新序列长度变为2N-1:

|n|≤N-1

设抗混叠插值滤波器为h(n)。

则零值内插后序列为:

v(n2Δ)=f(n/2Δ),n=0,±2,±4....0,otherwise

抗混叠滤波后得到:

f(n2Δ)=Σkv(k2Δ)h(n-k)=Σkf(kΔ)h(n-2k),

其中-(N-1)≤n≤N-1,-(N-1)/2≤k≤(N-1)/2。

r0(n)=h(2n+1)

r1(n)=h(2n)

当n=2m时:

f1(mΔ)=f(2m2Δ)=Σkf(kΔ)h(2m-2k)=Σkf(kΔ)r1(m-k)

当n=2m+1时:

f0(mΔ)=f(2m+12Δ)=Σkf(kΔ)h(2m+1-2k)=Σkf(kΔ)r0(m-k)

chirpA(n2Δ)=exp(-jπ4(Δ)2n2tanφ2),|n|≤N-1

chirpB(n2Δ)=exp(jπ4(Δ)2n2cscφ),|n|≤2N-1

其中φ≡aπ/2,a为变换阶次。上面两个chirp序列的多相分量分别如下:

e0(nΔ)=chirpB(2n2Δ),e1(nΔ)=chirpB(2n+12Δ)

l0(nΔ)=chirpA(2n2Δ),l1(nΔ)=chirpA(2n+12Δ)

则令

g(n2Δ)=chirpA(n2Δ)*f(n2Δ),|n|≤N-1

Ga(m2Δ)=Σkg(k2Δ)chirpB(m2Δ-k2Δ),|m|≤N-1

根据(5)式有

{Faf}(m2Δ)=Aφ2Δ*chirpA(m2Δ)*Ga(m2Δ),|m|≤N-1    (6)

把各个多相分量代入(6)式化简,并按两倍速率抽取{Faf}(m/(2Δ))即可得到输入信号序列f(n/Δ)的阶次为a的N点分数阶傅立叶变换结果{Faf}(m/Δ):

{Faf}(mΔ)=Aφ2Δ*l0(mΔ)*{Σk[Σkf(kΔ)r1(k-k)]l0(kΔ)e0(m-kΔ)

+Σk[Σkf(kΔ)r0(k-k)]l1(kΔ)e1(m-kΔ)},|m|≤(N-1)/2    (7)

其中|m|≤(N-1)/2,Aφ≡exp(-iπsgn(sinφ)/4+iφ/2)/|sinφ|。(7)式即可按前面介绍的5个步骤实现整个计算过程,整个过程具有两个支路并行运算的特点。

需要说明的是,在上述推导过程中并没有限定插值滤波器的类型,即本推导过程对任意一种插值滤波器都适用。

本发明提出的一种基于多抽样的分数阶傅立叶变换高效实现方法,其有益效果在于:

(1)本发明提出的实现方法相比于原始算法经典实现方法运算量更小。设插值滤波器长度为4N,那么两种实现方法的运算量如表1所示,从表1中可以看出多抽样高效实现方法每一支路的运算量均小于原始算法运算量的一半;

(2)本发明提出的实现方法在结构上具有两路并行的特点,实现时可以采用并行结构来提高计算速度,节省计算时间,尤其适合于FPGA等硬件实现。具体可参见“具体实施方式”中的FPGA实现方法实例;

(3)本发明提出的实现方法选用低阶次的加窗半带滤波器作为插值滤波器,节省了插值滤波环节一半的计算量,并且加窗可以有效防止在变换阶次较低时出现的振荡现象,具体见附图2中的比较。从图中可以看出插值滤波器在阶次较低,不加窗情况下,分数阶变换结果出现了振荡,误差大,而加窗后,则消除了振荡。

表1 本发明提出的实现方法和经典实现方法运算量的比较

  运算量  原始算法经典实现  多抽样实现单一支路  复乘数  12N log2N+36N  (11Nlog2N+22N)/2  复加数  24N log2N+48N  (22Nlog2N+23N)/2

附图说明

图1为多抽样高效实现方法流程图;

图2为Hanning窗半带滤波器与其它插值滤波器对分数阶傅立叶变换结果影响的比较;

图3是一个输入信号长度为1023点的本发明方法的FPGA实现结构框图;

具体实施方式

下面结合附图及FPGA实施例对发明内容做详细说明。

本发明涉及到一种基于多抽样的分数阶傅立叶变换高效实现方法,其原理见式(7),实现方法的算法流程图如图1所示,整个流程分解成以下5个步骤完成运算:

①将输入信号序列f(n/Δ)通过插值滤波器h(n)的多相分量ri(k)进行滤波,得到序列fi(m/Δ),其中i=0,1表示支路数。|n|≤(N-1)/2,Δ=N表示输入信号序列归一化的采样间隔,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);

②用chirp信号chirpA(n/(2Δ))的两个多相分量li(k/Δ)分别与步骤①相应支路的滤波结果fi(m/Δ)相乘得到gi(k/Δ),其中

chirpA(n2Δ)=exp(-jπ4(Δ)2n2tanφ2),|n|≤N-1,φ≡aπ/2,a为变换阶次;

l0(nΔ)=chirpA(2n2Δ),l1(nΔ)=chirpA(2n+12Δ);

③gi(k/Δ)与另一chirp信号chirpB(n/(2Δ))的两个多相分量ei(k/Δ)分别作卷积,其中

chirpB(n2Δ)=exp(jπ4(Δ)2n2cscφ),|n|≤2N-1;

e0(nΔ)=chirpB(2n2Δ),e1(nΔ)=chirpB(2n+12Δ);

④将步骤③中的两路卷积结果分别截取N点;

⑤将步骤④得到的两个支路的N点卷积结果相加,相加结果再与步骤②中的序列l0(k/Δ)和系数Aφ/(2Δ)相乘输出最后结果{Faf}(m/Δ),其中,|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滤波器就可实现,其中

h(n)=w(n)·2sin(πn2)πn,w(n)=0.5·[1+cos(L)],n=-L,…-1,0,1,…L,L=61

(2)ChirpDDS模块采用传统的Chirp DDS设计方法,由相位累加单元和频率累加单元再加上控制单元和正弦查找表构成。主要功能是输出公式(6)中对应的几个chirp序列:chirpA(n/2Δ),chirpB(n/2Δ)。序列chirpA的两个多相分量l0(k/Δ),l1(k/Δ)用于步骤②两个支路的调制。序列chirpB的两个多相分量e0(k/Δ),e1(k/Δ)用于步骤③的卷积运算,其中

chirpA(n2Δ)=exp(-jπ4(Δ)2n2tanφ2),|n|≤N-1;

chirpB(n2Δ)=exp(jπ4(Δ)2n2cscφ),|n|≤2N-1;

e0(nΔ)=chirpB(2n2Δ),e1(nΔ)=chirpB(2n+12Δ);

l0(nΔ)=chirpA(2n2Δ),l1(nΔ)=chirpA(2n+12Δ);

(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 xc2vp20FPGA芯片中实现,工作在100M时钟,那么整个计算过程消耗的时间大约为4.1ms,这要比传统方法的FPGA实现的9.5ms快一倍左右。

以上所述的具体描述,对发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限定本发明的保护范围,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号