CN1188938A - 一种用于谱分析中离散傅立叶变换的高精度数据处理方法 - Google Patents

Abstract

一种用于谱分析中离散傅立叶变换的高精度数据处理方法,属信号处理领域。该方法是在原连续Fourier变换的基础上,对离散后的样本函数进行精确积分,从而得到具有高精度的离散Fourier变换表达式,定义出“谱修正乘子”,在原FFT算法的基础上,只要乘上“谱修正乘子”就可得到高精度结果;本方法的结果无混迭效应、无Nyquist抽样率的限制,计算精度高、计算量不大。

Description

一种用于谱分析中离散傅立叶变换的高精度数据处理方法

本发明涉及一种用于谱分析的高精度数据处理方法，属于信号处理技术领域。

在当今的许多科学研究中，傅立叶(Fourier)变换实质上是一个分析和求解问题的普遍方法和重要工具，从基本性质看，Fourier变换的重要性还在于它使人们可以用一个全然不同的新观点去研究一些传统的领域。因此，在解决问题时，同时把一个函数及其Fourier变换确定出来，往往是求解问题的关键。

现有的用于谱分析的技术主要为经典的离散Fourier变换DFT(DiscreteFourier Transform)及其快速算法一快速Fourier变换FFT(Fast FourierTransform)，由于它的简单性和通用性，作为一种重要的分析手段，已广泛应用于线性系统、变换理论、信号处理、仿真、通信理论、光学、随机过程、概率论、量子物理等领域。在经典的离散Fourier变换DFT中，由于对原连续信号函数采取的是脉冲序列抽样函数并进行加窗处理，所得到的脉冲状的离散抽样函数和原始连续信号函数差别很大，因而所得到的Fourier变换存在精度差(特别是在较高频率的区域)及混迭效应现象。Nyquist抽样率就是为防止混迭效应现象所定义的最低抽样频率。因而，在实际的谱分析中，抽样频率一般都必须取得较高，即抽样点较多，因而FFT所进行的计算量很大，即使这样所得到的谱分析结果精度还是较差(只有较低阶的结果才满足要求)。

本发明之目的在于：提供一种用于谱分析中离散Fourier变换的高精度数据处理方法，基于该数据处理方法，可以设计出用于谱分析的专用仪器仪表、电子线路、专用芯片及硬件装置，或采用基于该数据处理方法编制的通用软件来实现。

本发明用于谱分析的离散Fourier变换的高精度数据处理方法，包括下述步骤：

(1)对连续信号函数f(t)进行离散化抽样，即 $\hat{f}\left(t\right)=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t\right)\mathrm{\δ}(t-\mathrm{k\Δt})$ 这里Δt为抽样间隔，δ(t)为脉冲函数。

在本发明中，每一抽样间隔内是用分段线性函数(折线)来逼近原连续信号函数f(t)。(2)对离散的抽样信号f(t_k)，首先按经典的FFT计算变换。 $\hat{F}\left({\mathrm{\ω}}_j\right)=\frac{\mathrm{\Δt}}{T_0}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_k}$ $=\frac{1}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_k}$ 以及首尾两点变换 $\hat{F}(t_0,t_N,{\mathrm{\ω}}_j)=\frac{\text{\Δt}}{T_0}[f\left(t_N\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_N}-f\left(t_0\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_0}]$ $=\frac{1}{N}[f\left(t_N\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_N}-f\left(t_0\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_0}]$ 其中T0为总抽样区间，i²＝-1

(3)对应于每一ω_j，定义相应的“谱修正乘子”C₁(ω_j)及C₂(ω_j)为 $C_1\left({\mathrm{\ω}}_j\right)={\left(\frac{\sin \left(\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}\right)}{\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\right)}^2$ $C_2\left({\mathrm{\ω}}_j\right)=\frac{1}{2}{\left(\frac{\sin \left(\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}\right)}{\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\right)}^2+\frac{1}{2\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\[1-\frac{\sin \left(2\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}\right)}{2\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\]i$ (4)将“谱修正乘子”C₁(ω_j)及C₂(ω_j)分别作用在经典FFT变换的结果

及首尾两点变换

上就可以得到具有高精度的离散Fourier变换，即 $F\left({\mathrm{\ω}}_j\right)=C_1\left({\mathrm{\ω}}_j\right)\·\hat{F}\left({\mathrm{\ω}}_j\right)+C_2\left({\mathrm{\ω}}_j\right)\·\hat{F}(t_0,t_N,{\mathrm{\ω}}_j)$ 该式可以计算对应于任意ω_j的Fourier变换。(5)该方法的主要计算量是在于 的计算，实际上这就是经典的Fourier变换，这可以用经典的FFT进行处理，对于具有N个抽样点的f(t_k)(k＝0，1，.......，N-1)，用FFT处理需要Nlog₂N次乘法， 的计算为2N次乘法，

的总计算为是3N次乘法，最后F(ω_j)的计算约为Nlog₂N+N+3N＝N(log₂N+4)次乘法。

(6)对应于高精度离散Fourier变换F(ω_j)的逆变换，可用Fourier级数的谐波合成定理来获取，为 $f\left(t_k\right)=\underset{j=-M}{\overset{M}{\mathrm{\Σ}}}F\left({\mathrm{\ω}}_j\right)\·e^{2\mathrm{\πi}{\mathrm{\ω}}_j{t}_k}$ M为谐波的阶数，M的选取视对逆变换的精度要求而定，可以是一个很大的数。

以上就连续函数给出了谱分析的新型高精度离散Fourier变换及快速算法，下面就一组给定的实或复的序列给出高精度离散Fourier变换方法及快速算法。

这时，可在以上连续Fourier变换的基础上取

t₀＝0 ${\mathrm{\ω}}_j=\frac{j}{N\•\mathrm{\Δt}}$        j＝0，1，…t_k＝k·Δt         k＝0，1，…N-1T₀＝N·Δt则对离散的序列f(t_k)(k＝1，2，3，….N-1)，经典的FFT公式为 $\hat{F}\left(j\right)=\frac{1}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)e^{-2\mathrm{\πijk}/N}$ 以及首尾两点变换 $\hat{F}(t_0,t_N,j)=\frac{\text{1}}{N}[f\left(t_N\right)-f\left(t_0\right)]$ “谱修正乘子”C₁(j)和C₂(j)分别为 $C_1\left(j\right)={\left(\frac{\sin (\mathrm{\π}j/N)}{\mathrm{\πj}/N}\right)}^2$ $C_2\left(j\right)=\frac{1}{2}{C}_1\left(j\right)+\frac{1}{2\mathrm{\πj}/N}\[1-\frac{\sin (2\mathrm{\πj}/N)}{2\mathrm{\πj}/N}\]i$ 那么，具有高精度的离散Fourier变换为 $F\left(j\right)=C_1\left(j\right)\·\hat{F}\left(j\right)+C_2\left(j\right)\·\hat{F}(t_0,t_N,j)$ j＝0，±1，±2，…该表达式对于任意的j都成立。对应于具有高精度的离散Fourier变换的逆变换为： $f\left(t_k\right)=\underset{j=-M}{\overset{M}{\mathrm{\Σ}}}F\left(j\right)\·e^{2\mathrm{\πijk}/N}$ M的选取视逆变换的精度要求而定，可以是一个很大的数。

同样，对于离散序列f(t_k)(k＝1，2，3，.......N-1)的高精度Fourier变换，如果

的计算采用FFT，其总计算量约为N(log₂N+4)次乘法。

下面给出逆变换的具体表达式：假设

         t₀＝0

       f(t_N)＝f(t₀)则 $f\left(t_k\right)=\underset{j=-M}{\overset{M}{\mathrm{\Σ}}}F\left(j\right)e^{2\mathrm{\πijk}/N}$ $=F\left(0\right)+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\left\{\right[F\left(j\right)+F(-j)]\cos (2\mathrm{\πjk}/N)-i[F\left(j\right)-F(-j)\left]\sin (2\mathrm{\πjk}/N)\right\}$ $=F\left(0\right)+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\left\{2C_1\left(j\right)\right[\frac{1}{N}\underset{l=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_l\right)\cos (2\mathrm{\πjl}/N)\left]\right\}\cos (2\mathrm{\πjk}/N)$ $+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\left\{2C_1\left(j\right)\right[\frac{1}{N}\underset{l=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_l\right)\sin (2\mathrm{\πjl}/N)\left]\right\}\sin (2\mathrm{\πjk}/N)$ $=F\left(0\right)+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}[a_j\cos (2\mathrm{\πjk}/N)+b_j\sin (2\mathrm{\πjk}/N)]$ $=F\left(0\right)+\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\left[A_j\sin (2\mathrm{\πjk}/N+{\mathrm{\φ}}_j)\right]$ $=F\left(0\right)+I_m\left[\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{A}_j{e}^{i(2\mathrm{\πjk}/N+{\mathrm{\φ}}_j)}\right]$ 其中： $a_j=R_e\left[2F\left(j\right)\right]=2C_1\left(j\right)\·R_e\left[\hat{F}\left(j\right)\right]$ $b_j=-I_m\left[2F\left(j\right)\right]=-2C_1\left(j\right)\·I_m\left[\hat{F}\left(j\right)\right]$ $A_j=\sqrt{a_j^2+b_j^2}$ ${\mathrm{\φ}}_j=\arctan \frac{a_j}{b_j}$

本发明的方法可用于各种数理及工程问题的谱分析，为一通用的高精度数据处理方法，采用该方法处理所得到的结果具有很高的分析精度，对真实的频谱具有很好的逼近性，并且无混迭效应、无抽样率的限制，可得到任意阶的谱分析结果，总体计算量不大。可根据实际情况选取大间隔进行抽样，这样可以大大减少计算量，特别是当原信号函数为分段折线时的情形，该谱分析方法可以得到精确的结果。

附图说明：

图1表示本发明的实施流程框图。

图2表示一连续函数f(t)＝e^-t(t＞0)的抽样点。

图3表示对连续函数f(t)＝e^-t(t＞0)进行离散Fourier变换的实部。

图4表示对连续函数f(t)＝e^-t(t＞0)进行离散Fourier变换的虚部。

图5表示通过对函数f(t)＝e^-t(t＞0)的离散Fourier变换来计算Fourier逆变换。

图6表示一组实序列f(t)(k＝0，1，2，......，20)。

图7表示对实序列f(t)进行离散Fourier变换的实部。

图8表示对实序列f(t)进行离散Fourier变换的虚部。

图9表示通过对实序列f(t)的离散Fourier变换来计算Fourier逆变换。

图10表示一矩形波的抽样点。

图11表示对矩形波采用经典Fourier变换和逆变换的结果。

图12表示对矩形波采用本发明方法进行Fourier变换和逆变换的结果。

结合附图说明实施例如下：

1、连续函数Fourier变换的离散计算

图2所示函数为f(t)＝e^-t(t＞0，当t＜0时此函数定义为零)，我们希望用离散方法来计算这个函数的Fourier变换的近似值。取 ${\mathrm{\ω}}_j=\frac{j}{N\·\mathrm{\Δt}}$ j＝0，1，…N-1t_k＝k·Δt该函数Fourier变换的解析结果为 $F\left(j\right)=\frac{1}{1+{(j/(N\·\mathrm{\Δt}))}^2}-\left[\frac{j/(N\·\mathrm{\Δt})}{1+{(j/(N\·\mathrm{\Δt}))}^2}\right]i$ 下面采用本发明的方法计算该函数的Fourier变换和逆变换：

(1)选择样本数N和抽样时间间隔Δt，图2中给出了函数f(t)＝e^-t(t＞0)在N＝40、Δt＝0.25时的样本，即

抽样A：N＝40、Δt＝0.25为使Fourier逆变换能够成立，在间断点上的函数值必须取中值，即取f(t＝0)≌0.5。

(2)对离散的抽样点，按经典的快速Fourier变换(FFT)计算变换 $\hat{F}\left(\frac{j}{N\·\mathrm{\Δt}}\right)=\mathrm{\Δt}\·\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}\left[e^{-k\·\mathrm{\Δt}}\right]\·e^{-2\mathrm{\πijk}/N}$ j＝0，1，2，…N-1注意，为了使连续和离散变换两者之间等价，在上式中，引入了比例因子Δt；按以上经典方法所得到的结果示于图3和图4中。

同时计算首尾两点变换 $\hat{F}(t_O,t_N,j)=\mathrm{\Δt}[f\left(t_N\right)e^{-2\mathrm{\πijk}/N}-f\left(t_O\right)e^{-2\mathrm{\πij}\·O/N}]$ ＝Δt[f(t_N)-f(t₀)](3)对应于各阶频率j，分别计算“谱修正乘子”C₁(j)和C₂(j)，即 $C_1\left(j\right)={\left(\frac{\sin \left(\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}\right)}{\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\right)}^2={\left(\frac{\sin (\mathrm{\πj}/N)}{\mathrm{\πj}/N}\right)}^2$ $C_2\left(j\right)=\frac{1}{2}{C}_1\left(j\right)+\frac{1}{2\mathrm{\πj}/N}[1-\frac{\sin (2\mathrm{\πj}/N)}{2\mathrm{\πj}/N}]i$ j＝0，1，2，…(4)将“谱修正乘子”C₁(j)和C₂(j)分别作用在经典FFT变换的结果

及首尾两点变换 上，即 $F\left(j\right)=C_1\left(j\right)\·\hat{F}\left(j\right)+C_2\left(j\right)\·\hat{F}(t_O,t_N,j)$ j＝0，1，2，…

表1列出了经典FFT方法、本发明的方法和解析结果的比较。图3和图4也表示了以上结果的比较。从结果可以看出，经典FFT方法的结果在高频区其误差很大，而本发明的方法无论在高频区还是在低频区都可以非常精确地逼近解析结果。

表示各种方法对连续函数f(t)＝e^-t(t＞0)进行离散

表1  Fourier变换的结果比较。

j	j/N	经典方法(FFT)	本发明方法(HFFT)	精确值(EXACT)
					02468101214161820222426284550556090100	.000.050.100.150.200.250.300.350.400.450.500.550.600.650.7001.1251.2501.3751.5002.2502.500	(.1130E+00， .0000E+00i)(.5179E-01，-.4807E-01i)(.2670E-01，-.3304E-01i)(.1962E-01，-.2279E-01i)(.1687E-01，-.1646E-01i)(.1556E-01，-.1212E-01i)(.1486E-01  -.8869E-02i)(.1445E-01，-.6245E-02i)(.1422E-01，-.3992E-02i)(.1409E-01，-.1948E-02i)(.1405E-01，-.9329E-09i)(.1409E-01， .1948E-02i)(.1422E-01， .3992E-02i)(.1445E-01， .6245E-02i)(.1486E-01， .8869E-02i)(.2224E-01，-.2725E-01i)(.1556E-01，-.1212E-01i)(.1432E-01，-.5084E-02i)(.1405E-01， .1474E-09i)(.1556E-01，-.1212E-01i)(.1405E-01，-.7611E-08i)	(.1005E+00， .0000E+00i)(.3897E-01，-.4898E-01i)(.1374E-01，-.3453E-01i)(.6606E-02，-.2491E-01i)(.3825E-02，-.1924E-01i)(.2481E-02，-.1561E-01i)(.1736E-02，-.1311E-01i)(.1280E-02，-.1129E-01i)(.9827E-03，-.9907E-02i)(.7774E-03，-.8826E-02i)(.6300E-03，-.7957E-02i)(.5204E-03，-.7244E-02i)(.4368E-03，-.6649E-02i)(.3713E-03，-.6145E-02i)(.3188E-03，-.5713E-02i)(.1141E-03，-.3502E-02i)(.9926E-04，-.3171E-02i)(.8308E-04，-.2889E-02i)(.6999E-04，-.2652E-02i)(.3063E-04，-.1764E-02i)(.2520E-04，-.1591E-02i)	(.1000E+00， .0000E+00i)(.3877E-01，-.4872E-01i)(.1367E-01，-.3435E-01i)(.6574E-02，-.2478E-01i)(.3807E-02，-.1914E-01i)(.2470E-02，-.1552E-01i)(.1729E-02，-.1303E-01i)(.1276E-02，-.1122E-01i)(.9798E-03，-.9850E-02i)(.7757E-03，-.8773E-02i)(.6293E-03，-.7908E-02i)(.5206E-03，-.7197E-02i)(.4378E-03，-.6602E-02i)(.3733E-03，-.6098E-02i)(.3221E-03，-.5666E-02i)(.1249E-03，-.3532E-02i)(.1012E-03，-.3180E-02i)(.8367E-04，-.2891E-02i)(.7031E-04，-.2651E-02i)(.3126E-04，-.1768E-02i)(.2532E-04，-.1591E-02i)

(5)为说明采用本发明的方法无混迭效应、无Nyquist抽样率的限制，下面再分别采用另外二种抽样率对f(t)＝e^-t(t＞0)进行抽样，并分别计算相应的Fourier变换，即

抽样B：N＝80、Δt＝0.125

抽样C：N＝160、Δt＝0.0625几种抽样方案的有关结果及其比较列于表2中，结果表明：除鞍点(即对应于ω_j＝j/(N·Δt)＝1/Δt，或j＝N的频率)外，对应于任意阶ω_j的变换都非常精确地逼近解析结果，并且无混迭效应和Nyquist抽样率的限制。由于抽样点多，则抽样函数可以更好地逼近原函数，所以计算的精度自然要高，当然计算量也相应增加。

表示在三种抽样率下，采用本发明方法对连续函数

表2  f(t)＝e^-t(t＞0)进行离散Fourier变换的结果比较。

j	N＝40  Δt＝0.25	N＝80  Δt＝0.125	N＝160  Δt＝0.0625
				02468101214161820222426284550556090100	(.1005E+00， .0000E+00 i)(.3897E-01，-.4898E-01 i)(.1374E-01，-.3453E-01 i)(.6606E-02，-.2491E-01 i)(.3825E-02，-.1924E-01 i)(.2481E-02，-.1561E-01 i)(.1736E-02，-.1311E-01 i)(.1280E-02，-.1129E-01 i)(.9827E-03，-.9907E-02 i)(.7774E-03，-.8826E-02 i)(.6300E-03，-.7957E-02 i)(.5204E-03，-.7244E-02 i)(.4368E-03，-.6649E-02 i)(.3713E-03，-.6145E-02 i)(.3188E-03，-.5713E-02 i)(.1141E-03，-.3502E-02 i)(.9926E-04，-.3171E-02 i)(.8308E-04，-.2889E-02 i)(.6999E-04，-.2652E-02 i)(.3063E-04，-.1764E-02 i)(.2520E-04，-.1591E-02 i)	(.1001E+00， .0000E+00 i)(.3882E-01，-.4878E-01 i)(.1368E-01，-.3439E-01 i)(.6582E-02，-.2481E-01 i)(.3812E-02，-.1916E-01 i)(.2473E-02  -.1554E-01 i)(.1731E-02，-.1305E-01 i)(.1277E-02，-.1124E-01 i)(.9809E-03，-.9862E-02 i)(.7766E-03，-.8785E-02 i)(.6300E-03，-.7918E-02 i)(.5212E-03，-.7206E-02 i)(.4383E-03，-.6611E-02 i)(.3737E-03，-.6107E-02 i)(.3223E-03，-.5674E-02 i)(.1249E-03，-.3538E-02 i)(.1011E-03，-.3186E-02 i)(.8348E-04，-.2897E-02 i)(.6999E-04，-.2657E-02 i)(.3054E-04，-.1764E-02 i)(.2520E-04，-.1590E-02 i)	(.1000E+00， .0000E+00 i)(.3878E-01，-.4874E-01 i)(.1367E-01，-.3436E-01 i)(.6575E-02，-.2479E-01 i)(.3808E-02，-.1914E-01 i)(.2471E-02，-.1553E-01 i)(.1729E-02，-.1304E-01 i)(.1276E-02，-.1123E-01 i)(.9800E-03，-.9853E-02 i)(.7759E-03，-.8776E-02 i)(.6294E-03，-.7910E-02 i)(.5208E-03，-.7199E-02 i)(.4380E-03，-.6604E-02 i)(.3734E-03，-.6100E-02 i)(.3221E-03，-.5667E-02 i)(.1250E-03，-.3533E-02 i)(.1012E-03，-.3181E-02 i)(.8368E-04，-.2892E-02 i)(.7032E-04，-.2652E-02 i)(.3126E-04，-.1769E-02 i)(.2532E-04，-.1592E-02 i)

(6)通过对函数f(t)＝e^-t(t＞0)的高精度离散Fourier变换应用谐波合成定理来计算Fourier逆变换。结果见图5。2、一组序列的Fourier变换

对于一组实序列f(k)(k＝1，2，……，20)(见图6)，即

5，10，-8，3，-2，1，18，-3，1，-8，4，-5，-1，6，-10，10，-1，4，-2，3，5按以下步骤计算该序列的Fourier变换：

(1)按经典的快速Fourier变换(FFT)计算变换 $\hat{F}\left(j\right)=\frac{1}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(k\right)e^{-2\mathrm{\πijk}/N}$ j＝0，1，2，…N-1这时首尾两点变换 $\hat{F}(t_O,t_N,j)=0$ (2)对应于各阶j，分别计算“谱修正乘子”C₁(j)，即 $C_1\left(j\right)={\left(\frac{\sin (\mathrm{\πj}/N)}{\mathrm{\πj}/N}\right)}^2$ j＝0，1，2，…(3)将“谱修正乘子”C₁(j)作用在经典FFT变换的结果

上，即 $F\left(j\right)=C_1\left(j\right)\·\hat{F}\left(j\right)$

这就是本发明所提出的高精度Fourier变换方法的结果，表3列出了经典FFT方法和本发明方法对序列f(k)的变换结果及其比较。实际上，对形如图6所示的样本函数的Fourier变换，表3所列的结果是精确的，即：采用本发明的方法对分段折线的信号函数的处理是精确的。

表示各种方法对一组实序列f(k)进行离散Fourier变换的表3

结果比较。

j	j/N	经典方法(FFT)	本发明方法(HFFT)
				012345	.000.050.100.150.200.250	( .1250E+01，  .0000E+00 i)( .8295E+00， -.3737E+00 i)(-.6113E+00，  .5584E+00 i)( .1223E+01，  .5264E+00 i)( .1077E+01， -.1384E+01 i)(-.2390E-06， -.2500E+00 i)	( .1250E+01， .0000E+00 i)( .8227E+00，-.3706E+00 i)(-.5914E+00， .5403E+00 i)( .1135E+01， .4886E+00 i)( .9427E+00，-.1211E+01 i)(-.2811E-06，-.2026E+00 i)

67891011121314151617181920

.300.350.400.450.500.550.600.650.700.750.800.850.900.9501.000

( .7863E+00，   .9035E+00 i)( .7556E-01，  -.1441E+01 i)( .7977E+00，  -.1907E+01 i)(-.1878E+01，   .1959E+01 i)(-.8500E+00，  -.2272E-05 i)(-.1878E+01，  -.1959E+01 i)( .7977E+00，   .1907E+01 i)( .7556E-01，   .1441E+01 i)( .7863E+00，  -.9035E+00 i)( .1370E-07，   .2500E+00 i)( .1077E+01，   .1384E+01 i)( .1223E+01，  -.5264E+00 i)(-.6113E+00，  -.5584E+00 i)( .8295E+00，   .3737E+00 i)( .1250E+01，   .2371E-05 i)

( .5794E+00， .6657E+00 i)( .4961E-01，-.9462E+00 i)( .4569E+00，-.1092E+01 i)(-.9169E+00， .9562E+00 i)(-.3445E+00，-.7680E-06 i)(-.6138E+00，-.6401E+00 i)( .2031E+00， .4855E+00 i)( .1438E-01， .2743E+00 i)( .1064E+00，-.1223E+00 i)( .1657E-07， .2252E-01 i)( .5892E-01， .7572E-01 i)( .3536E-01，-.1522E-01 i)(-.7302E-02，-.6670E-02 i)( .2278E-02， .1027E-02 i)(-.2403E-06， .7220E-20 i)

(4)通过对实序列f(k)的高精度离散Fourier变换应用谐波合成定理来计算Fourier逆变换。即： $f_k=F\left(0\right)+I_m\left[\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{A}_j{e}^{i(2\mathrm{\πjk}/N+{\mathrm{\φ}}_j)}\right]$ 其中 $A_j=\sqrt{a_j^2+b_j^2}$ ${\mathrm{\φ}}_j=\arctan \frac{a_j}{b_j}$ $a_j=2c_1\left(j\right)\·R_e\text{[}\hat{F}\left(j\right)\text{]}$ $b_j=-2C_1\left(j\right)\·I_m\left[\hat{F}\left(j\right)\right]$ 结果见表4。表4

表示用实序列f(k)的离散Fourier变换来计算Fourier逆变换。

序列	本发明的逆变换(inverse HFFT)	实序列点
			12345678	( .5071E+01.     .0000E+00 i)( .9426E+01，    .0000E+00 i)(-.7276E+01，    .0000E+00 i)( .2604E+01，    .0000E+00 i)(-.1799E+01，    .0000E+00 i)( .1347E+01，    .0000E+00 i)( .1705E+02，    .0000E+00 i)(-.2380E-01，    .0000E-00 i)	.5000E+01.1000E+02-.8000E+01.3000E+01-.2000E+01.1000E+01.1800E+02-.3000E+0，

9101112131415161718192021

( .6808E+00，    .0000E+00 i)(-.7476E+01，    .0000E+00 i)( .3477E+01，    .0000E+00 i)(-.4676E+01，    .0000E+00 i)(-.9261E+00，    .0000E+00 i)( .5428E+01，    .0000E+00 i)(-.9103E+01，    .0000E+00 i)( .9228E+01，    .0000E+00 i)(-.6044E+00，    .0000E+00 i)( .3728E+01，    .0000E+00 i)(-.1725E+01，    .0000E+00 i)( .2927E+01，    .0000E+00 i)( .5071E+01，    .0000E+00 i)

.1000E+01-.8000E+01.4000E+01-.5000E+01-.1000E+01.6000E+01-.1000E+02.1000E+02-.1000E+01.4000E+01-.2000E+01.3000E+01.5000E+01

图10、图11和图12表示对另一矩形波采用经典方法和本发明方法进行Fourier变换和逆变换的结果。

Claims (1)

1.一种用于谱分析中离散Fourier变换的高精度数据处理方法，其特征在于包括下述步骤：

(1)对连续信号函数f(t)进行离散化抽样，即 $\hat{f}\left(t\right)=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t\right)\mathrm{\δ}(t-\mathrm{k\Δt})$ 这里Δt为抽样间隔，δ(t)为脉冲函数。

在本发明中，每一抽样间隔内是用分段线性函数(折线)来逼近原连续信号函数f(t)。

(2)对离散的抽样信号f(t_k)，首先按经典的FFT计算变换。 $\hat{F}\left({\mathrm{\ω}}_j\right)=\frac{\mathrm{\Δt}}{T_0}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_k}$ $=\frac{1}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_k}$ 以及首尾两点变换 $\hat{F}(t_0,t_N,{\mathrm{\ω}}_j)=\frac{\text{\Δt}}{T_0}[f\left(t_N\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_N}-f\left(t_0\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_0}]$ $=\frac{1}{N}[f\left(t_N\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_N}-f\left(t_0\right)e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_0}]$ (3)对应于每一ω_j，定义相应的“谱修正乘子”C₁(ω_j)及C₂(ω_j)为 $C_1\left({\mathrm{\ω}}_j\right)={\left(\frac{\sin \left(\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}\right)}{\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\right)}^2$ $C_2\left({\mathrm{\ω}}_j\right)=\frac{1}{2}{\left(\frac{\sin \left(\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}\right)}{\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\right)}^2+\frac{1}{2\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\[1-\frac{\sin \left(2\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}\right)}{2\mathrm{\π}{\mathrm{\ω}}_j\mathrm{\Δt}}\]i$ (4)将“谱修正乘子”C₁(ω_j)及C₂(ω_j)分别作用在经典FFT变换的结果

及首尾两点变换 上就可以得到具有高精度的离散Fourier变换，即 $F\left({\mathrm{\ω}}_j\right)=C_1\left({\mathrm{\ω}}_j\right)\·\hat{F}\left({\mathrm{\ω}}_j\right)+C_2\left({\mathrm{\ω}}_j\right)\·\hat{F}(t_0,t_N,{\mathrm{\ω}}_j)$ 该式可以计算对应于任意ω_j的Fourier变换。

(5)对应于高精度离散Fourier变换F(ω_j)的逆变换，可用Fourier级数的谐波合成定理来获取，为 $f\left(t_k\right)=\underset{j=-M}{\overset{M}{\mathrm{\Σ}}}F\left({\mathrm{\ω}}_j\right)\·e^{-2\mathrm{\πi}{\mathrm{\ω}}_j{t}_k}$ M的选取视对逆变换的精度要求而定，可以是一个很大的数。

(6)对于一组给定的实或复的序列f(t_k)(k＝1，2，3，.…N-1)，它的高精度离散Fourier变换方法及快速算法为：

经典的FFT公式为 $\hat{F}\left(j\right)=\frac{1}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)e^{-2\mathrm{\πijk}/N}$ 以及首尾两点变换 $\hat{F}(t_0,t_N,j)=\frac{1}{N}[f\left(t_N\right)-f\left(t_0\right)]$ “谱修正乘子”C₁(j)和C₂(j)分别为 $C_1\left(j\right)={\left(\frac{\sin (\mathrm{\πj}/N)}{\mathrm{\πj}/N}\right)}^2$ $C_2\left(j\right)=\frac{1}{2}{C}_1\left(j\right)+\frac{1}{2\mathrm{\πj}/N}\[1-\frac{\sin (2\mathrm{\πj}/N)}{2\mathrm{\πj}/N}\]i$ 那么，具有高精度的离散Fourier变换为 $F\left(j\right)=C_1\left(j\right)\·\hat{F}\left(j\right)+C_2\left(j\right)\·\hat{F}(t_0,t_N,j)$ j＝0，±1，±2，…该表达式对于任意的j都成立。

对应于具有高精度的离散Fourier变换的逆变换为： $f\left(t_k\right)=\underset{j=-M}{\overset{M}{\mathrm{\Σ}}}F\left(j\right)\·e^{2\mathrm{\πijk}/N}$ M的选取视逆变换的精度要求而定，可以是一个很大的数。

同样，对于离散序列f(t_k)(k＝1，2，3……N-1)的高精度Fourier变换，逆变换的具体表达式可表达为： $f\left(t_k\right)=\underset{j=-M}{\overset{M}{\mathrm{\Σ}}}F\left(j\right)e^{2\mathrm{\πijk}/N}$ $=F\left(0\right)+I_m\left[\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{A}_j{e}^{i(2\mathrm{\πjk}/N+{\mathrm{\φ}}_j)}\right]$ 其中： $a_j=R_e\left[2F\left(j\right)\right]=2C_1\left(j\right)\·R_e\left[\hat{F}\left(j\right)\right]$ $b_j=-I_m\left[2F\left(j\right)\right]=-2C_1\left(j\right)\·I_m\left[\hat{F}\left(j\right)\right]$ $A_j=\sqrt{a_j^2+b_j^2}$ ${\mathrm{\φ}}_j=\arctan \frac{a_j}{b_j}$

Images