fft窄带高分辨率算法

fft窄带高分辨率算法,第1张

fft算法是什么

FFT算法(fast Fourier transform),即快速傅里叶变换,是指利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称,简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。

概念

有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列。但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT)。 1965年,Cooley和Tukey提出了计算离散傅里叶变换(DFT)的快速算法,将DFT的运算量减少了几个数量级。从此,对快速傅里叶变换(FFT)算法的研究便不断深入,数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT的多种算法,基本算法是基2DIT和基2DIF。FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。

快速傅氏变换(FFT),是离散傅氏变换的快速算法,它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换,可以说是进了一大步。

设x(n)为N项的复数序列,由DFT变换,任一X(m)的计算都需要N次复数乘法和N-1次复数加法,而一次复数乘法等于四次实数乘法和两次实数加法,一次复数加法等于两次实数加法,即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法),那么求出N项复数序列的X(m),即N点DFT变换大约就需要N^2次运算。当N=1024点甚至更多的时候,需要N2=1048576次运算,在FFT中,利用WN的周期性和对称性,把一个N项序列(设N=2k,k为正整数),分为两个N/2项的子序列,每个N/2点DFT变换需要(N/2)^2次运算,再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后,总的运算次数就变成N+2*(N/2)^2=N+N^2/2。继续上面的例子,N=1024时,总的运算次数就变成了525312次,节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去,直到分成两两一组的DFT运算单元,那么N点的DFT变换就只需要Nlog2N次的运算,N在1024点时,运算量仅有10240次,是先前的直接算法的1%,点数越多,运算量的节约就越大,这就是FFT的优越性。

如何提高fft算法分辨率

FFT程序,输入是一组复数,输出也是一组复数,想问一下输入到底应该输入什么,输出的复数的含义是什么。

给定一组序列的抽样值,如何用FFT确定它的频率。

首先,fft函数出来的应该是个复数,每一个点分实部虚部两部分。

假设采用1024点fft,采样频率是fs,那么第一个点对应0频率点,第512点对应的就是fs/2的频率点。然后从头开始找模值最大的那个点,其所对应的频率值应该就是要的基波频率了。

FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。

有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。

这就是很多信号分析采用FFT变换的原因。

另外,FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。

虽然很多人都知道FFT是什么,可以用来做什么,怎么去做,但是却不知道FFT之后的结果是什么意思、如何决定要使用多少点来做FFT。

一个模拟信号,经过ADC采样之后,就变成了数字信号。

采样定理告诉,采样频率要大于信号频率的两倍,这些就不在此罗嗦了。

采样得到的数字信号,就可以做FFT变换了。

N个采样点,经过FFT之后,就可以得到N个点的FFT结果。

为了方便进行FFT运算,通常N取2的整数次方。

在地球物理数据处理中,经常遇到处理二维实数据的情况。例如在地震勘探中,对面波勘探数据作频散分析解释时,要将时间-空间域的信息转换为频率-波数域频谱;在重磁异常的滤波或转换中,要将空间域的异常f(x,y)转换为波数域F(ω,υ)等。这些分析都需要进行二维的傅里叶变换(FFT)。

根据傅里叶变换的定义,对于连续二维函数f(x,y),其傅里叶变换对为

地球物理数据处理基础

对于离散的二维序列fjk(j=0,1,…,M-1;k=0,1,…,N-1),其傅里叶变换为

地球物理数据处理基础

1.二维复序列的FFT算法

对于M条测线,每条测线N个测点,构成复序列yjk(j=0,1,…,M-1;k=0,1,…,N-1),根据离散傅里叶公式(8-41),其傅里叶变换为

地球物理数据处理基础

于是,可以分两步套用一维复FFT完成二维复FFT的计算。

(1)沿测线方向计算

对于j=0,1,…,M-1逐测线套用一维复FFT,执行式(8-43)。定义复数组 则算法为

1)对于j=0,1,…,M-1,作2)~7);

2)将yjk输入A1(k),即A1(k)=yjk(k=0,1,…,N-1);

3)计算Wr,存入W(r),即

4)q=1,2,…,p(p=log2N),若q为偶数执行6),否则执行5);

5)k=0,1,2,…,(2p-q-1)和n=0,1,2,…,(2q-1-1)循环,作

A2(k2q+n)=A1(k2q-1+n)+A1(k2q-1+n+2p-1)

A2(k2q+n+2q-1)=[A1(k2q-1+n)-A1(k2q-1+n+2p-1)]·W(k2q-1)

至k,n循环结束;

6)k=0,1,2,…,(2p-q-1)和n=0,1,2,…,(2q-1-1)循环,作

A1(k2q+n)=A2(k2q-1+n)+A2(k2q-1+n+2p-1)

A1(k2q+n+2q-1)=[A2(k2q-1+n)-A2(k2q-1+n+2p-1)]·W(k2q-1)

至k,n循环结束;

7)q循环结束,若p为偶数,将A1(n)输入到Yjn,否则将A2(n)输入到Yjn(n=0,1,…,N-1);

8)j循环结束,得到Yjn(j=0,1,…,M-1;n=0,1,…,N-1)。

(2)垂直测线方向计算

对于n=0,1,…,N-1逐一套用一维复FFT,执行式(8-44)。即

1)对于n=0,1,…,N-1,作2)~7);

2)将Yjn输入A1(j),即A1(j)=Yjn(j=0,1,…,M-1);

3)计算Wr存入W(r),即

4)q=1,2,…,p(p=log2M),若q为偶数执行6),否则执行5);

5)j=0,1,2,…,(2p-q-1)和m=0,1,2,…,(2q-1-1)循环,作

A2(j2q+m)=A1(j2q-1+m)+A1(j2q-1+m+2p-1)

A2(j2q+m+2q-1)=[A1(j2q-1+m)-A1(j2q-1+m+2p-1)]·W(j2q-1)

至j,m循环结束;

6)j=0,1,2,…,(2p-q-1)和m=0,1,2,…,(2q-1-1)循环,作

A1(j2q+m)=A2(j2q-1+m)+A2(j2q-1+m+2p-1)

A1(j2q+m+2q-1)=[A2(j2q-1+m)-A2(j2q-1+m+2p-1)]·W(j2q-1)

至j,m循环结束;

7)q循环结束,若p为偶数,将A1(m)输入到Ymn,否则将A2(m)输入到Ymn(m=0,1,…,M-1);

8)n循环结束,得到二维复序列的傅氏变换Ymn(m=0,1,…,M-1;n=0,1,…,N-1),

所求得的Ymn是复数值,可以写为

Ymn=Rmn+iImn (m=0,1,…,M-1;n=0,1,…,N-1)

其中,Rmn,Imn的值也是已知的。

2.二维实序列的FFT算法

对于二维的实序列,我们把其看作是虚部为零的复序列,套用上述的二维复序列FFT方法来求其频谱算法上也是可行的,但势必会增加大量的无功运算。因此,有必要研究二维实序列FFT的实用算法,同一维实序列FFT的实现思路一样,同样把二维实序列按一定的规律构造成二维复序列,调用二维复序列FFT,然后通过分离和加工得到原实序列的频谱。

例如采样区域有2 M条测线,每条测线有N个点,并且M,N都是2的整数幂,需要计算实样本序列xjk(j=0,1,2,…,2 M-1;k=0,1,2,…,N-1)的傅氏变换:

地球物理数据处理基础

类似于一维实序列FFT的思想,直接建立下面的二维实序列FFT算法:

(1)将一个二维实序列按偶、奇线号分为两个二维子实序列,分别作为实部和虚部组合为一个二维复序列。即令

地球物理数据处理基础

(2)调用二维复FFT过程,求出yjk的二维傅氏变换Ymn的复数值:

地球物理数据处理基础

式中:Rmn,Imn是Ymn的实部和虚部。

(3)利用Rmn,Imn换算Xmn的值。

前两步容易实现,下面分析第(3)步的实现。

记hjk,gjk的傅氏变换为Hmn,Gmn。根据傅里叶变换的定义,我们导出Xmn与Hmn,Gmn的关系式:

地球物理数据处理基础

式中,Hmn,Gmn为复数,我们用上标r和i表示其实部和虚部,将上式右端实部、虚部分离

地球物理数据处理基础

其中:

地球物理数据处理基础

下面的任务是将Hmn,Gmn各分量与通过二维复FFT求出的Rmn,Imn值联系起来。为此先给出奇、偶分解性质和类似于一维情况的三个二维傅氏变换性质:

(1)奇偶分解性

任何一个正负对称区间定义的函数,均可唯一地分解为如下偶(even)、奇(odd)函数之和:

地球物理数据处理基础

(2)周期性

地球物理数据处理基础

(3)复共轭性

地球物理数据处理基础

现在我们来建立Rmn,Imn与Hmn,Gmn的关系。对Ymn作奇偶分解:

地球物理数据处理基础

根据线性性质

地球物理数据处理基础

对照式(8-54)和式(8-55),得

地球物理数据处理基础

由于hjk,gjk是实函数,根据复共轭性质,上面两式对应的奇偶函数相等。即

地球物理数据处理基础

再由奇偶分解性和周期性,得

地球物理数据处理基础

将式(8-57)代入式(8-50),得

地球物理数据处理基础

再利用Hmn,Gmn周期性及复共轭性,可以得到m=M/2+1,…,M-1;n=0,1,…,N-1的傅氏变换,即

地球物理数据处理基础

将式(8-50)中M,N改为M-m,N-n,并将上式代入,得

地球物理数据处理基础

由式(8-58)、式(8-59)和式(8-61)即可得到原始序列xjk(j=0,1,…,2M-1;n=0,1,…,N-1)在m=0,1,…,M-1;n=0,1,…,N-1区间的傅氏变换Xmn。

具体二维实序列的FFT算法如下:

(1)令hjk=x2j,k,gjk=x2j+1,k,形成

yjk=hjk+igjk (j=0,1,…,2 M-1;n=0,1,…,N-1)

(2)调用二维复序列FFT过程,即从两个方向先后调用一维复FFT算法式(8-43)和式(8-44),求得yjk的二维傅氏变换Ymn的复数值:

Ymn=Rmn+iImn (m=0,1,…,M-1;n=0,1,…,N-1)

(3)用下列公式由Rmn,Imn的值换算Xmn的值:

地球物理数据处理基础

地球物理数据处理基础


欢迎分享,转载请注明来源:内存溢出

原文地址:https://www.54852.com/yw/7612992.html

(0)
打赏 微信扫一扫微信扫一扫 支付宝扫一扫支付宝扫一扫
上一篇 2023-04-07
下一篇2023-04-07

发表评论

登录后才能评论

评论列表(0条)

    保存