CN105137176B - 一种利用快速三角形式傅里叶变换的信号谐波分析方法 - Google Patents
一种利用快速三角形式傅里叶变换的信号谐波分析方法 Download PDFInfo
- Publication number
- CN105137176B CN105137176B CN201510493899.XA CN201510493899A CN105137176B CN 105137176 B CN105137176 B CN 105137176B CN 201510493899 A CN201510493899 A CN 201510493899A CN 105137176 B CN105137176 B CN 105137176B
- Authority
- CN
- China
- Prior art keywords
- base
- harmonic
- time
- decimation
- signal
- 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.)
- Active
Links
Abstract
本发明公开一种利用快速三角形式傅里叶变换的信号谐波分析方法,包括如下步骤:对信号进行采样,得到N=2L点信号序列;对信号序列进行基2i时间抽取,i为不大于L‑1的正整数,得到2i个2L‑i点信号序列;对各2L‑i点信号序列进行傅里叶变换,得到各2L‑i点序列的自身所有的谐波系数;将各2L‑i点序列的各次谐波系数变换成N点信号序列从1到2L‑i‑1次谐波过程数据;将基从2i变换至基2i‑1,并不断地重复变换过程直至基为2,并计算出对应于N点信号序列的各次谐波系数,由此得到信号的傅里叶表达式。本发明所述技术方案显著提高了对信号进行傅里叶变换的计算效率,提高了信号谐波分析的实时性能。
Description
技术领域
本发明涉及信号谐波的分析方法。更具体地,涉及一种利用快速三角形式傅里叶变换(Fast Triangular-form Fourier Transform,FTFT)的信号谐波分析方法。
背景技术
傅里叶变换在物理学、电子类学科、数论、组合数学、信号处理、概率论、统计学、密码学、声学、光学、海洋学、结构动力学等领域都有着广泛的应用。然而,由于傅里叶变换计算复杂,很大程度上制约了傅里叶的应用。1965年,Cooley和Tukey在《计算机科学》发表著名的《机器计算傅立叶级数的一种算法》论文,快速傅里叶变换(Fast Fourier Transform-FFT)开始大规模应用。Cooley和Tukey的FFT算法的最基本运算为蝶形运算,每个蝶形运算包括两个输入点,因而也称为基2算法。在这之后,又有一些新的算法,进一步提高了FFT的运算效率,比如基4算法,分裂基算法等,这些新算法对FFT运算效率的提高一般在50%以内,因此,效率提高的效果有限。上述所有计算方法都是在傅里叶变换的复指数形式上进行的,其计算推理过程相对复杂,也不利于理解。
因此,需要提供一种快速的傅里叶变换方法应用于数字信号处理,以期提高傅里叶变换的效率,实现快速对信号的谐波进行分析,提升数字设备的信号实时处理性能。
发明内容
本发明的目的在于提供一种利用快速三角形式傅里叶变换(Fast Triangular-form Fourier Transform,FTFT)的信号谐波分析方法,可以快速获得信号各次谐波的系数。在系数计算过程中,无须引入复数的计算形式,变换过程容易理解。采用发明中提出变基时间抽取法(谐波衍生法),无需进行FFT运算中对数据进行的倒位序排列,即可快速计算各次谐波的系数,较大的提高了计算效率,使乘法计算量最大减少到基2FFT的25%以下。另外,利用变基进行谐波分析的方法具有较大的灵活性,可以据此衍生出各种定基、混合基等计算方法,使信号谐波分析更加方便。
为达到上述目的,本发明采用下述技术方案:
一种利用快速三角形式傅里叶变换的信号谐波分析方法,该方法包括如下步骤:
对信号进行周期性采样,每周期得到N点信号序列,并满足N=2L;
对N点信号序列进行基2i时间抽取,i为小于L的正整数,得到2i个2L-i点信号序列;
对各2L-i点信号序列进行傅里叶变换,得到各2L-i点序列的自身所有的谐波系数;
将各2L-i点序列的各次谐波系数变换成N点信号序列从1到2L-i-1次谐波过程数据;
利用变基时间抽取的方法特性,将基从2i变换至基2i-1,并不断地重复变换过程直至基为2,并计算出对应于N点信号序列的各次谐波系数,由此得到信号的傅里叶表达式。
优选地,步骤“将各2L-i点序列的各次谐波系数变换成N点信号序列从1到2L-i-1次谐波过程数据”的计算方法:
式中,Rnij表示第n谐波基2i时间抽取时的第j个余弦系数的过程量;Inij为第n谐波基2i时间抽取时的第j个正弦系数的过程量;R'nij表示基2i时间抽取后,第j个子序列第n次谐波余弦系数;I'nij表示j个子序列第n次谐波的正弦系数;α=2π/N,N为T周期内信号的采样点数;jα为时间抽取后第j个子序列第n次谐波系数向合成序列的第j个过程量转换的转换基角,其中j的取值为区间[0,2i-1)的整数。
优选地,步骤“利用变基时间抽取的方法特性,将基从2i变换至基2i-1,并不断地重复变换过程直至基为2,并计算出对应于N点信号序列的各次谐波系数,由此得到信号的傅里叶表达式”中计算基2i-1时间抽取对应的各次谐波系数的计算方法为:
用于频率衍生的计算方法:
用于迭代过程的计算方法:
式中,Rnij表示第n谐波基2i时间抽取时的第j个余弦系数的过程量;Inij为第n谐波基2i时间抽取时的第j个正弦系数的过程量;Yij为基2i时间抽取后,子序列所有数据的代数和。
优选地,利用变基时间抽取法实现基2、4或8谐波抽取的信号谐波分析方法,该方法包括如下步骤:
对信号进行采样,得到N点信号序列,按基2Q,Q取1、2或3抽取的方法,获取最小序列点数为2Q的子序列,并计算各子序列的所有谐波系数;
将顺序相邻的2Q个子序列组通过变基计算获取包含22Q个数据的子序列的所有谐波系数,此时子序列个数变为2L-2Q,转换基角为2πj/22Q。
不断的重复变基计算合成过程,使合成序列包含的数据点数不断增加,直到合成的数列点数为N。
优选地,当采样时间点数N不能使log4 N和log8 N为整数时,利用改变时间抽取的最小序列点数,实现基4、8时间抽取的倒位序方法:
对应于原序列序号为m的数据,用(bL-1bL-2…b1b0)表示m值,则:
对应于最小序列点数为2的基8倒位序方法(b2b1b0b5b4b3…bL-2bL-3bL-4bL-1);
对应于最小序列点数为4的基8倒位序方法(b2b1b0b5b4b3…bL-3bL-4bL-5bL-1bL-2);
对应于最小序列点数为2的基4倒位序方法(b1b0b3b2…bL-2bL-3bL-1);
上述计算方法中采用基8时,用于第一次计算的最小子序列的点数不为8,用于基4计算的最小子序列点数为2。
优选地,当采样时间点数N不能使log4 N和log8 N为整数时,采用混合基方法实现信号谐波分析的倒位序方法,方法实现基8基4混合,基8基2混合,基4基2混合的方法:
对应于原序列序号为m的数据,用(bL-1bL-2…b1b0)表示m值,则:
基8基4混合的倒位序方法为(b1b0b4b3b2…bL-1bL-2bL-3);
基8基2混合的倒位序方法为(b0b3b2b1…bL-1bL-2bL-3);
基4基2混合的倒位序方法为(b0b2b1…bL-1bL-2);
上述方法中基8计算方法的最后一次的计算中采用的基不为8,而基4方法的除最后一次计算采用的基不为4。
一种利用快速三角形式傅里叶变换的信号谐波分析方法中的谐波变换方法,在实现基从2i变换至基2i-1,部分实现部分谐波从基2i-1至2i-2的变换,该方法的变换公式为:
式中,表示第N/2i+n次谐波基2i-1时间抽取时第j个余弦系数过程量,为第N/2i+n次谐波基2i-1时间抽取时第j个正弦系数过程量。
本发明的有益效果如下:
本发明所述技术方案在计算过程中,无须引入复数的计算形式,变换过程容易理解;采用谐波系数推衍法,快速计算各次谐波的系数,较大的提高了计算效率,提高了信号处理的实时性能。例如,在对1024点信号序列进行傅里叶变换时,利用变基FTFT方法的乘法计算次数可降至采用基2FFT算法的所需乘法计算次数25%以下。在信号周期内的采样数据点数越小,相对提升的效果就越显著。
附图说明
下面结合附图对本发明的具体实施方式作进一步详细的说明。
图1示出基于变基时间抽取的利用快速三角形式傅里叶变换的信号谐波分析方法流程图。
图2示出基于变基时间抽取的利用快速三角形式傅里叶变换的信号谐波分析方法中谐波系数的中间过程数据的获取示意图。
图3示出基于定基时间抽取的利用快速三角形式傅里叶变换的信号谐波分析方法流程图。
图4示出基于混合基时间抽取的利用快速三角形式傅里叶变换的信号谐波分析方法流程图。
具体实施方式
为了更清楚地说明本发明,下面结合优选实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解,下面所具体描述的内容是说明性的而非限制性的,不应以此限制本发明的保护范围。
对于任意连续信号f(t),其三角形式的傅里叶级数可以用式(1)表示:
式(1)中:
式(2)中,T为信号f(t)的周期,t0表示计时起点,n表示第n次谐波,ω=2π/T。可将式(2)离散化即可得到对应傅里叶级数的离散形式。在T周期内,对信号进行N次采样,获得对应于信号在T周期内的N点序列,式(3)给出了信号序列的标识形式。
XN={x0,x1,x2,...,xN-1} (3)
由于N点序列最多可计算至N/2次谐波的系数,因此对于N点序列,其对应的傅里叶级数可以表示为:
式(4)中:
在式(5)中α0=ωt0,α=ωT/N=2π/N,且:
基于变基时间抽取的利用快速三角形式傅里叶变换的信号谐波分析方法,包括如下步骤:
对信号进行周期性采样,每周期得到N点信号序列,并满足N=2L,L为正整数;
对N点信号序列进行基2i时间抽取,i为小于L的正整数,得到2i个2L-i点信号序列;
对各2L-i点信号序列进行傅里叶变换,得到各2L-i点序列的自身所有的谐波系数;
将各2L-i点序列的各次谐波系数变换成N点信号序列从1到2L-i-1次谐波过程数据;
利用变基时间抽取的方法特性,将基从2i变换至基2i-1,并不断地重复变换过程直至基为2,并最后计算出对应于N点信号序列的各次谐波系数,由此得到信号的傅里叶表达式。
该方法的原理为:
式(5)给出了一般离散傅里叶变化的三角函数形式。由此可知,要计算各次谐波的系数值,只需计算Rn与In的值,即各次谐波的余弦和正弦系数。为了简化Rn和In的计算,可对数列进行基2i时间抽取,此时N点序列被分成2i个N/2i点序列,对Rn计算可得:
同理可推得In的表达式,因此可获得如式所示的傅里叶计算式:
式(8)中:
式(8)表明,可以将Rn及In的计算分成2i个段,并对各个段的结果进行求和即可得到N点序列第n次谐波的系数值。将式(9)展开可得:
在式(10)中:
对比式(6)可知,式(11)实际上是以2i为基对序列进行时间抽取,获取2i个数列后,对应于第j个序列傅里叶变换各次谐波系数的计算式,而式(10)实际上给出了如何将第j个序列的第n次谐波系数转换为N点序列第n次谐波系数第j个段的计算结果表达式,jα为第j个序列基波系数的转换角度,第j个序列其余各次谐波,其系数转换角度均为jα的整数倍,因此,可以将这个角度定义为对应于第j个序列的转换基角。由α的值可知,对N点序列进行基2i抽取后,对应于第j个序列的谐波系数转换基角为2πj/N。由此可知,转换基角与抽取前序列的点数N相关。
式(9)~(11)说明Rn的计算值可从抽取后对应数据系列的谐波系数获得。在进行基2i时间抽取后,原序列的数据个数(点数)N减小为N/2i,可以起到简化计算的目的。另一方面,由于抽取后序列的数据个数为N/2i,因此,每个序列实际上最多可计算的谐波个数为N/2i+1。当谐波次数大于N/2i+1时,必须通过其他方法获得。由式(6)可计算N/2i-n的谐波系数表达式,其中计算过程如式(12)所示
同理可推得的计算式,并通过与式(11)对比可知式(13)成立:
式(13)实际上给出了N/2i-n次谐波在基2i-1抽取时,各段数据与n次谐波系数在基2i时间抽取时系数的关系,在采用基2i时间抽取的方法对序列进行整理计算后,只需获得前N/2i+1次谐波得中间过程数据,并保留谐波计算过程中对应的计算过程量,即可从这些过程量中计算N/2i-n次谐波的系数。
另外,由式(11)可以推得基2i时间抽取法的中间过程量与基2i-1抽取计算过程量的关系,其中Rn(i-1)j的计算过程如下:
同理可求得In(i-1)j与Inij之间的关系,因此,式(14)成立
由式(14)可知,在获得基2i时间抽取后,要获得基2i-1时间抽取的前N/2i+1次只需将间隔为2i-1的谐波系数对应两两相加即可。式(13)和式(14)中,j为区间[0,2i-1)内的整数,进行式(13)和式(14)的整理后,可以获得从1到N/2i-1次谐波系数的过程量。
进行上述整理后,各个谐波的j值减小了一半,然而谐波次数增加了约一倍,因此整个计算过程中所保存的中间过程量基本不变。这个过程相当于从基2i时间抽取变换成基2i-1时间抽取的过程,因此,可以称为变基计算过程。然而上述变换过程不是完整的,因为基2i-1时间抽取时,必须获得N/2i次谐波的计算过程量,而由式(13)和式(14)无法获得。此时,可将n=N/2i直接代入式(9),并应用式(14)所述的关系可得:
式(15)中
由式(15)可知Yij为对N点序列进行基2i谐波抽取后,对应于第j个序列包含的所有数据和,因此,在计算过程中应保存抽取后所有各个序列所有数据的代数和,从而方便计算能被2i整除的特殊次谐波。另外,由式(16)可知:
式(17)说明了n=N/2i与n=N/2i-1谐波系数的中间过程关系。这两个系数具有递推的关系,即当获得n=N/2i对应的Y值时,要计算n=N/2i-1的谐波可以从n=N/2i的谐波过程中间结果继续,这与式(14)的计算推演过程类似。
式(13)~(15)给出了完整的从基2i到基2i-1的变换过程,在这个过程中,参与计算的谐波个数增加了一倍,而各次谐波系数的中间过程量减少了一倍,因此,总的中间过程数据量不变。
由上述分析可知,变基时间抽取的计算方法可以用图1表示:
由图1可知变基时间主要可以归结为三级骤,其中第1级完成对N点序列基2i的时间抽取,在第1级计算中,求出抽取后各子序列的谐波系数是计算量所在。当抽取后子序列的点数为4时,子序列谐波系数的计算不包含乘法计算,从而可以使第一级计算简化。第二级完成后将各子序列的谐波系数转化为N点序列前2L-i-1次谐波过程数据,第三级计算完成基变换,并计算出最后谐波系数的结果。
下面以初始基为2L-2为例,说明该方法的具体实现过程:
第一级计算:
对于任意N(N>4)点序列,当满足N=2L时,总可以对其进行基2L-2时间抽取,此时,N点序列被分成2L-2个4点序列。
计算各个4点序列的谐波系数,其中1次谐波系数包含R和I两个值,而2次谐波系数仅包含R,I为0,另外,必须保存所有4点序列对应的数据和,用于下次迭代过程。
第二级计算:
利用式(10)整理出基2L-2时间抽取的1、2次谐波的计算中间过程量,并令x=2
第三级计算:
利用式(13)计算2x-1+1到2x-1次谐波系数基2L-x-1中间过程量,其中n取所有1到2x -1-1之间的值
利用式(15)计算2x次基2L-x-1谐波系数中间过程量。
利用式(14)计算出N点序列1到2x-1次谐波系数基2L-x-1的中间过程量,并用(17)整理2x+1次谐波需要的Y值。
x值加1,若x值小于L则持续第三级计算过程,继续进行变基计算,否则结束计算过程。
在上述计算过程中,x与上述推导过程中i的关系为i=L-x,此处选用x主要是为了计算过程中方便表示。
在上面的示例中,最大时间基(即第一次计算的时间基)为2L-2,由于基于时间抽取后每个序列谐波系数计算过程都是完整的,这里选用2L-2是因为利用这种方式抽取后,所有抽取后的序列仅包含4点数据,其谐波系数计算不包含乘法运算,较有代表性。
在上述方法中第三级谐波的计算中,谐波系数的中间过程数据的获取可以图2来表达。
从图2可以看出,大部分增加的谐波过程数据都可以从上一次谐波的基础上衍生出来,例如7次谐波从1次谐波衍生,6次谐波从二次谐波数据衍生等等。从这个意义上说,变基方法也可以成为谐波衍生方法。式(13)实际上构成了谐波衍生的依据,即在原有谐波的基础上,通过一定的方法得出相关的谐波的过程数据,并用式(15)增加一项无法从现有谐波数据衍生的谐波。从式(13)可知,若原有谐波的基为2i,则衍生谐波的基为2i-1,即衍生谐波的基总比原谐波的基小一倍。此时,原有谐波再通过式(14)将基变为2i-1,并应用式(16)为下一次增加的谐波做数据处理,实现变基过程。在这个过程中,原谐波的过程数据减少一半,但谐波的数量增加了一倍,因此整个过程中,用于计算的过程数据总量是不变的。
以下为该方法的计算量分析:
在基于时间抽取的傅里叶变换计算过程中,第1级不包含乘法运算,第二级计算中,第一个序列由于变换角度为0,无需进行乘法计算,因此参与式(10)计算的序列共有2L -2-1个。在计算1次谐波时,每个序列1次谐波系数需进行4次乘法计算,因此,1次谐波系数转换共需4×(2L-2-1)次乘法计算,而计算2次谐波时,由于I值为0,因此,包含2×(2L-2-1)次乘法计算,由此可知,第二级计算共包含6(2L-2-1)次乘法计算。
在第三级计算中,每个谐波的j值需取遍0到2L-x-1-1的所有整数值,共2L-x-1个数,由于式(13)中每个j值需进行4次乘法计算,而当j取0和2L-x-2,对应的角度值为0和π/2,不包含乘法计算,当j取2L-x-3和3×2L-x-3时,对应的角度值为π/4和3π/4,可简化为2次乘法计算,此时,相当于j减少3个取值的乘法运算,而对于式(15)当j取0和2L-x-2时不包含乘法计算,其余数值只包含2个乘法计算。由于每个x值有2x-1-1个谐波需进行式(13)计算以及1个谐波需进行式(15)计算,因此,对于任意x值,第4、5步总的乘法计算量为:
4×(2x-1-1)(2L-x-1-3)+2×(2L-x-1-2)=2L-3×2x+1-2L-x+8
由于x取L-1和L-2时(即满足2L-x-1-3>0),第三级运算过程不包含乘法计算,由此可知上述方法总的乘法计算次数为:
需要注意的是,在上面的计算过程中,计算次数第二项存在的条件是L-3≥2,即L≥5,这说明只有32点及以上的序列才会有第二项计算次数,同时也说明8、16点序列采用这种方法时具有相对较高的效率。实际上通过合理的优化,对于8点序列最终只需要2次乘法计算,而16点序列可以简化至12次乘法计算。这也说明如果采用的基为2L-3或者2L-4等较大的时间基将会有更高的计算效率,这里不再举例说明。
由上述乘法计算公式可知,对于1024点的序列,进行上述变基过程之后,所需总的乘法次数仅为5706,由于基2FFT计算方法的复数乘法计算次数约为(N/2)log2 N,而一次复数乘法包含4次实数乘法,因此其实数乘法的计算次数约为2N log2 N,对于1024点序列,基2FFT约需20480次实数乘法。因此,上述变基时间抽取方式的乘法运算仅为基2FFT的27.86%。
另外,上述算法的效率仍可以继续提高。式(13)给出了N/2i-n与n次谐波的计算关系,为从2i时间基向2i-1的时间基转换过程奠定了基础。类似的,可以推得N/2i+n与n次谐波的计算关系,如式(18)所示
对比式(13)可知,式(18)所需与式(13)所需进行的4个乘法计算项完全相同,因此计算了N/2i-n的谐波系数后,无需再对N/2i+n的谐波系数进行乘法计算。即运用式(13)实现N/2i-n次谐波基2i-1时间中间计算结果时,实际上也包含了一部分N/2i+n,基2i-1时间中间计算结果。例如,从1次谐波衍生3次谐波的数据时,5次谐波的过程数据实际上无需进行乘法计算。由此可知,从2i-1时间基向2i-2时间基转换的过程中,需要进行乘法运算的谐波个数减少了。如果用Mx表示x次计算过程需计算的谐波个数,则第x次与x-1次谐波计算的个数关系为:
Mx=2x-1-Mx-1 (19)
利用迭代法将式(19)进行展开可得:
当x为偶数和奇数时,Mx的计算结果为:
由于第x次原来需要进行的乘法次数为2x-1,对比式(20)可知,利用上述关系大约可减小第三级计算过程约1/3的计算量,对于1024点的序列,利用上述关系可以使乘法计算次数进一步减少到5000次以下,使计算效率得到更大的提高。
关于定基时间抽取的谐波系数计算方法:
在FFT中一般都采用定基时间抽取的方法,而采用的基一般为8,4,2。在定基计算时需要对原始数列进行倒位序排列,例如,对应于原始序列序号为m的数据,若用二进制数(bL-1bL-2…b1b0)表示m值,则在基2时间抽取时,其对应的序号为(b0b1…bL-2bL-1),基4时间抽取时其序号为(b1b0b3b2…bL-1bL-2),而基8时间抽取时为其序号为(b2b1b0b5b4b3…bL-1bL- 2bL-3)。需要说明的是基4和基8时间抽取法分别要求L能被2和3整除。
上述倒位序的过程也就是按基时间抽取的过程,其实质是,不断的对N点序列,以及抽取后的子序列进行基2、4、8时间抽取。变基计算方法也完全适用于这种模式。
变基方法的实质是将基于时间抽取的子序列谐波数据合成原序列谐波数据的过程,即通过式(10)将子序列的谐波系数整理成对应于子序列的合成序列的中间过程数据,进而由变基计算取得合成数列谐波系数。如图3所示,如果将子序列再次进行时间抽取,并且使抽取的规则相同时,即构成了定基时间抽取的过程。此时,子序列谐波系数计算的过程也可以进行变基计算。例如对于信号采样点数为512数据序列,进行3级时间抽取后,子序列的点数为8。首先直接计算各子序列的谐波系数(也可由变基时间抽取法计算),从而获得各子序列4次谐波系数,而后,顺序取每8个子序列为一组,由式(10)合成64点序列的前4次谐波系数过程值,此时,需注意式(10)的转换基角为2πj/64。通过变基计算方法,在获得所有8组序列的全部32次谐波系数值后,继续通过式(10)合成512点的前32次谐波系数过程值,此时,式(10)的转换基角为2πj/512,进而继续通过变基计算,得出所有谐波系数。
从上面可以看出,对于基取2Q(Q取1,2,3)的定基时间抽取法,其计算可以用用下面的过程表述:
首先将信号采样获得的N点序列,按基2Q抽取的方法,获取最小序列点数为2Q的子序列,并计算各子序列的所有谐波系数。
通过图1中第二级以及第三级完成基2Q变基计算,将顺序相邻的2Q个子序列合成的含22Q个数据序列的所有谐波系数。由于此时合成序列的点数为22Q,因此,第二级计算在应用式(10)时,转换基角为2πj/22Q。继续应用变基方法,利用顺序相邻的2Q个含22Q个数据的子序列计算含23Q个数据序列的所有谐波系数,变基计算中,式(10)的转换基角变为2πj/23Q。此后不断的重复变基计算过程,直到合成的数列点数为N。
在定基时间抽取方法中,基8抽取法的计算效率最高,基4抽取发次之,基2抽取法最慢。当采用基4、8抽取时,要求信号周期内采样的点数必须满足或者这对采样的点数构成了较大的限制,为了避免这种限制,可通过最小子序列的点数进行调节。其对应的倒位序规则续作小量变动,对应于原序列序号为m的数据,若用(bL-1bL-2…b1b0)表示m值,则:
在基8算法中若信号的采样点数N仅满足N=2L,当L被3除余2时,可采用的倒位序排列为(b2b1b0b5b4b3…bL-3bL-4bL-5bL-1bL-2),即最小序列点数为4,当被3除余1时,倒位序排列为(b2b1b0b5b4b3…bL-2bL-3bL-4bL-1),即最小序列点数为2。在进行第二级计算时,需注意每个序列的转换基角仅与合成后序列的点数有关。
同理在基4时间抽取时,L被2除余1时,可以采用的倒位序规则为(b1b0b3b2…bL- 2bL-3bL-1),即最小的序列点数为2。
此外,如图4所示也可以采用混合基计算方式,如对于基8计算方法中,当L被3除余2或者1时,可采用的倒位序规则分别为:
(b1b0b4b3b2…bL-1bL-2bL-3)和(b0b3b2b1…bL-1bL-2bL-3)
此时,最后1次子序列合成计算中,分别采用基4和基2的变基时间抽取法,由于最后1次子序列合成计算所采用的基与开始不同,因此称为混合基计算方法。
同理在基4时间抽取时,当L被2除余1时,可采用的倒位序规则为:
(b0b2b1…bL-1bL-2)
在最后一步合成计算时,采用基2方法计算。
显然,本发明的上述实施例仅仅是为清楚地说明本发明所作的举例,而并非是对本发明的实施方式的限定,对于所属领域的普通技术人员来说,在上述说明的基础上还可以做出其它不同形式的变化或变动,这里无法对所有的实施方式予以穷举,凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。
Claims (7)
1.一种利用快速三角形式傅里叶变换的信号谐波分析方法,其特征在于,该方法包括如下步骤:
对信号进行周期性采样,每周期得到N点信号序列,并满足N=2L;
对N点信号序列进行基2i时间抽取,i为小于L的正整数,得到2i个2L-i点信号序列;
对各2L-i点信号序列进行傅里叶变换,得到各2L-i点序列的自身所有的谐波系数;
将各2L-i点序列的各次谐波系数变换成N点信号序列从1到2L-i-1次谐波过程数据;
利用变基时间抽取的方法特性,将基从2i变换至基2i-1,并不断地重复变换过程直至基为2,并计算出对应于N点信号序列的各次谐波系数,由此得到信号的傅里叶表达式。
2.根据权利要求1所述的利用快速三角形式傅里叶变换的信号谐波分析方法,其特征在于,步骤“将各2L-i点序列的各次谐波系数变换成N点信号序列从1到2L-i-1次谐波过程数据”的计算方法:
式中,Rnij表示第n谐波基2i时间抽取时的第j个余弦系数的过程量;Inij为第n谐波基2i时间抽取时的第j个正弦系数的过程量;R'nij表示基2i时间抽取后,第j个子序列第n次谐波余弦系数;I'nij表示j个子序列第n次谐波的正弦系数;α=2π/N,N为T周期内信号的采样点数;jα为时间抽取后第j个子序列第n次谐波系数向合成序列的第j个过程量转换的转换基角,其中j的取值为区间的整数。
3.根据权利要求1所述的利用快速三角形式傅里叶变换的信号谐波分析方法,其特征在于,步骤“利用变基时间抽取的方法特性,将基从2i变换至基2i-1,并不断地重复变换过程直至基为2,并计算出对应于N点信号序列的各次谐波系数,由此得到信号的傅里叶表达式”中计算基2i-1时间抽取对应的各次谐波系数的计算方法为:
用于频率衍生的计算方法:
用于迭代过程的计算方法:
式中,Rnij表示第n谐波基2i时间抽取时的第j个余弦系数的过程量,表示第(N/2i-n)谐波基2i-1时间抽取时的第j个余弦系数的过程量,表示第n谐波基2i时间抽取时的第(j+2i-1)个余弦系数的过程量,表示第(N/2i)谐波基2i-1时间抽取时的第j个余弦系数的过程量,Rn(i-1)j表示第n谐波基2i-1时间抽取时的第j个余弦系数的过程量;Inij为第n谐波基2i时间抽取时的第j个正弦系数的过程量,为第(N/2i-n)谐波基2i-1时间抽取时的第j个正弦系数的过程量,为第n谐波基2i时间抽取时的第(j+2i-1)个正弦系数的过程量,为第(N/2i)谐波基2i-1时间抽取时的第j个正弦系数的过程量,In(i-1)j为第n谐波基2i-1时间抽取时的第j个正弦系数的过程量;Yij为基2i时间抽取后,对应于第j个序列包含的所有数据的代数和;为基2i时间抽取后,对应于第(j+2i -1)个序列包含的所有数据的代数和;Y(i-1)j为基2i-1时间抽取后,对应于第j个序列包含的所有数据的代数和。
4.根据权利要求1所述的利用快速三角形式傅里叶变换的信号谐波分析方法,其特征在于,利用变基时间抽取法实现基2、4或8谐波抽取的信号谐波分析方法,该方法包括如下步骤:
对信号进行采样,得到N点信号序列,按基2Q,Q取1、2或3抽取的方法,获取最小序列点数为2Q的子序列,并计算各子序列的所有谐波系数;
将顺序相邻的2Q个子序列组通过变基计算获取包含22Q个数据的子序列的所有谐波系数,此时子序列个数变为2L-2Q,转换基角为2πj/22Q;
不断的重复变基计算合成过程,使合成序列包含的数据点数不断增加,直到合成的数列点数为N。
5.根据权利要求4所述的利用快速三角形式傅里叶变换的信号谐波分析方法,其特征在于,当采样时间点数N不能使log4N和log8N为整数时,利用改变时间抽取的最小序列点数,实现基4、8时间抽取的倒位序方法:
对应于原序列序号为m的数据,用二进制数(bL-1bL-2…b1b0)表示m值,则:
对应于最小序列点数为2的基8倒位序方法(b2b1b0b5b4b3…bL-2bL-3bL-4bL-1);
对应于最小序列点数为4的基8倒位序方法(b2b1b0b5b4b3…bL-3bL-4bL-5bL-1bL-2);
对应于最小序列点数为2的基4倒位序方法(b1b0b3b2…bL-2bL-3bL-1);
上述计算方法中采用基8时,用于第一次计算的最小子序列的点数不为8,用于基4计算的最小子序列点数为2。
6.如权利要求4所述的利用快速三角形式傅里叶变换的信号谐波分析方法,其特征在于,当采样时间点数N不能使log4N和log8N为整数时,采用混合基方法实现信号谐波分析的倒位序方法,方法实现基8基4混合,基8基2混合,基4基2混合的方法:
对应于原序列序号为m的数据,用二进制数(bL-1bL-2…b1b0)表示m值,则:
基8基4混合的倒位序方法为(b1b0b4b3b2…bL-1bL-2bL-3);
基8基2混合的倒位序方法为(b0b3b2b1…bL-1bL-2bL-3);
基4基2混合的倒位序方法为(b0b2b1…bL-1bL-2);
上述方法中基8计算方法的最后一次的计算中采用的基不为8,而基4方法的除最后一次计算采用的基不为4。
7.一种利用快速三角形式傅里叶变换的信号谐波分析方法中的谐波变换方法,其特征在于,在实现基从2i变换至基2i-1,部分实现部分谐波从基2i-1至2i-2的变换,该方法的变换公式为:
式中,表示第N/2i+n次谐波基2i-1时间抽取时第j个余弦系数过程量,为第N/2i+n次谐波基2i-1时间抽取时第j个正弦系数过程量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510493899.XA CN105137176B (zh) | 2015-08-12 | 2015-08-12 | 一种利用快速三角形式傅里叶变换的信号谐波分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510493899.XA CN105137176B (zh) | 2015-08-12 | 2015-08-12 | 一种利用快速三角形式傅里叶变换的信号谐波分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105137176A CN105137176A (zh) | 2015-12-09 |
CN105137176B true CN105137176B (zh) | 2017-12-19 |
Family
ID=54722592
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510493899.XA Active CN105137176B (zh) | 2015-08-12 | 2015-08-12 | 一种利用快速三角形式傅里叶变换的信号谐波分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105137176B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106053940B (zh) * | 2016-08-09 | 2018-04-10 | 重庆大学 | 一种基于方波傅里叶级数分解的谐波分析方法 |
CN106374802B (zh) * | 2016-09-19 | 2019-01-25 | 上海新时达电气股份有限公司 | 电机驱动器的死区补偿电压值的自动调整方法 |
CN113567788A (zh) * | 2021-07-30 | 2021-10-29 | 北京易艾斯德科技有限公司 | 混合基fft在电力系统中的应用方法、装置、设备和介质 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN1062793A (zh) * | 1990-12-25 | 1992-07-15 | 能源部电力科学研究院 | 快速高精度信号频谱分析方法 |
CN101061474A (zh) * | 2004-06-10 | 2007-10-24 | 哈桑·塞希托格鲁 | 信号处理的矩阵定值方法和装置 |
CN101551790A (zh) * | 2008-04-03 | 2009-10-07 | 中兴通讯股份有限公司 | 快速傅立叶变换实现方法及装置 |
EP2113777A1 (en) * | 2008-05-02 | 2009-11-04 | Tektronix International Sales GmbH | Signal analyzer and method for displaying frequency domain data |
CN104699657A (zh) * | 2013-12-06 | 2015-06-10 | 中国科学院电子学研究所 | 一种快速实现傅里叶变换的方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP2005055412A (ja) * | 2003-08-01 | 2005-03-03 | High Speed Signal Processing Laboratory Inc | Fft型周波数解析装置。 |
WO2013173975A1 (zh) * | 2012-05-22 | 2013-11-28 | 深圳市英威腾电气股份有限公司 | 谐波检测方法及相关装置 |
-
2015
- 2015-08-12 CN CN201510493899.XA patent/CN105137176B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN1062793A (zh) * | 1990-12-25 | 1992-07-15 | 能源部电力科学研究院 | 快速高精度信号频谱分析方法 |
CN101061474A (zh) * | 2004-06-10 | 2007-10-24 | 哈桑·塞希托格鲁 | 信号处理的矩阵定值方法和装置 |
CN101551790A (zh) * | 2008-04-03 | 2009-10-07 | 中兴通讯股份有限公司 | 快速傅立叶变换实现方法及装置 |
EP2113777A1 (en) * | 2008-05-02 | 2009-11-04 | Tektronix International Sales GmbH | Signal analyzer and method for displaying frequency domain data |
CN104699657A (zh) * | 2013-12-06 | 2015-06-10 | 中国科学院电子学研究所 | 一种快速实现傅里叶变换的方法 |
Also Published As
Publication number | Publication date |
---|---|
CN105137176A (zh) | 2015-12-09 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
GB2547816B (en) | Actually-measured marine environment data assimilation method based on sequence recursive filtering three-dimensional variation | |
Gehrmann et al. | Two-loop master integrals for $ q\overline {q}\to VV $: the planar topologies | |
Aistleitner et al. | A central limit theorem for Latin hypercube sampling with dependence and application to exotic basket option pricing | |
CN105137176B (zh) | 一种利用快速三角形式傅里叶变换的信号谐波分析方法 | |
CN107241081B (zh) | 余弦调制滤波器组的稀疏fir原型滤波器的设计方法 | |
CN109117455A (zh) | 计算装置及方法 | |
CN103093085B (zh) | 基于典型相关分析的稳态诱发电位的分析方法 | |
Hu et al. | A novel generic fast Fourier transform pruning technique and complexity analysis | |
CN113657937A (zh) | 基于eemd-cnn+sae-rfr混合算法的日前电价预测方法 | |
Zhang et al. | Load prediction based on hybrid model of VMD-mRMR-BPNN-LSSVM | |
Bao et al. | LSFQ: A low precision full integer quantization for high-performance FPGA-based CNN acceleration | |
CN106982045A (zh) | 一种基于socp优化的eir‑cmfb结构的设计方法 | |
CN115829157A (zh) | 基于变分模态分解和Autoformer模型的化工水质指标预测方法 | |
CN104820581B (zh) | 一种fft和ifft逆序数表的并行处理方法 | |
Patel et al. | Time frequency analysis: a sparse S transform approach | |
WO2018170400A1 (en) | Apparatus and methods of providing an efficient radix-r fast fourier transform | |
CN109857980A (zh) | 一种快速傅里叶分析算法 | |
CN103870437A (zh) | 数字信号处理装置及其处理方法 | |
Heng-Li et al. | INTEGRATING EMD, CHAOS-BASED NEURAL NETWORK AND PSO FOR FINANCIAL TIME SERIES FORECASTING. | |
Okazaki et al. | A one dimensional mapping method for time series data | |
Malathi et al. | Review on fast complex multiplication algorithms and implementation | |
CN116881665A (zh) | 一种基于CMOA优化的TimesNet-BiLSTM光伏出力预测方法 | |
CN105278923A (zh) | 一种基于aic信息准则的信号源个数估计硬件电路及其实现方法 | |
AKANBI et al. | A RESIDUE NUMBER SYSTEM BASED CONVOLUTION NEURAL NETWORK ALGORITHM | |
CN115829135A (zh) | 基于注意力机制的燃气负荷预测方法、装置、设备及介质 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |