fft频谱算法
⑴ 二维实序列的快速傅里叶变换(FFT)
在地球物理数据处理中,经常遇到处理二维实数据的情况。例如在地震勘探中,对面波勘探数据作频散分析解释时,要将时间-空间域的信息转换为频率-波数域频谱;在重磁异常的滤波或转换中,要将空间域的异常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的值:
地球物理数据处理基础
地球物理数据处理基础
⑵ FFT原理的FFT基本原理
FFT是一种DFT的高效算法,称为快速傅立叶变换(fast Fourier transform)。FFT算法可分为按时间抽取算法和按频率抽取算法,先简要介绍FFT的基本原理。从DFT运算开始,说明FFT的基本原理。
DFT的运算为:
式中
由这种方法计算DFT对于X(K)的每个K值,需要进行4N次实数相乘和(4N-2)次相加,对于N个k值,共需N*N乘和N(4N-2)次实数相加。改进DFT算法,减小它的运算量,利用DFT中
的周期性和对称性,使整个DFT的计算变成一系列迭代运算,可大幅度提高运算过程和运算量,这就是FFT的基本思想。
FFT基本上可分为两类,时间抽取法和频率抽取法,而一般的时间抽取法和频率抽取法只能处理长度N=2^M的情况,另外还有组合数基四FFT来处理一般长度的FFT 设N点序列x(n),,将x(n)按奇偶分组,公式如下图
改写为:
一个N点DFT分解为两个 N/2点的DFT,继续分解,迭代下去,其运算量约为
其算法有如下规律
两个4点组成的8点DFT
四个2点组成的8点DFT
按时间抽取的8点DFT
原位计算
当数据输入到存储器中以后,每一级运算的结果仍然储存在同一组存储器中,直到最后输出,中间无需其它存储器
序数重排
对按时间抽取FFT的原位运算结构,当运算完毕时,这种结构存储单元A(1)、A(2),…,A(8)中正好顺序存放着X(0),X(1),X(2),…,X(7),因此可直接按顺序输出,但这种原位运算的输入x(n)却不能按这种自然顺序存入存储单元中,而是按X(0),X(4),X(2),X(6),…,X(7)的顺序存入存储单元,这种顺序看起来相当杂乱,然而它也是有规律的。当用二进制表示这个顺序时,它正好是“码位倒置”的顺序。
蝶形类型随迭代次数成倍增加
每次迭代的蝶形类型比上一次蝶代增加一倍,数据点间隔也增大一倍 频率抽取2FFT算法是按频率进行抽取的算法。
设N=2^M,将x(n)按前后两部分进行分解,
按K的奇偶分为两组,即
得到两个N/2 点的DFT运算。如此分解,并迭代,总的计算量和时间抽取(DIT)基2FFT算法相同。
算法规律如下:
蝶形结构和时间抽取不一样但是蝶形个数一样,同样具有原位计算规律,其迭代次数成倍减小 时,可采取补零使其成为
,或者先分解为两个p,q的序列,其中p*q=N,然后进行计算。 前面介绍,采用FFT算法可以很快算出全部N点DFT值,即z变换X(z)在z平面单位圆上的全部等间隔取样值。实际中也许①不需要计算整个单位圆上z变换的取样,如对于窄带信号,只需要对信号所在的一段频带进行分析,这时希望频谱的采样集中在这一频带内,以获得较高的分辨率,而频带以外的部分可不考虑,②或者对其它围线上的z变换取样感兴趣,例如语音信号处理中,需要知道z变换的极点所在频率,如极点位置离单位圆较远,则其单位圆上的频谱就很平滑,这时很难从中识别出极点所在的频率,如果采样不是沿单位圆而是沿一条接近这些极点的弧线进行,则在极点所在频率上的频谱将出现明显的尖峰,由此可较准确地测定极点频率。③或者要求能有效地计算当N是素数时序列的DFT,因此提高DFT计算的灵活性非常有意义。
螺旋线采样是一种适合于这种需要的变换,且可以采用FFT来快速计算,这种变换也称作Chirp-z变换。
⑶ 如何采用fft分析信号的频谱
你的问题太抽象了,如果是周期信号的话,一般是由现实的某个系统产生的,你可以估计它的周期,比如测G值的扭秤信号,大概一小时左右,那么你就可以根据抽样定理对它进行2倍于最大频率的抽样,再用fft进行频谱分析。如果是平稳随机信号,如语音信号你就知道它的频率主要分量在20-3k之内,你就可以用大于6k的对它采样,再进行fft。如果是非平稳信号,你就不能用fft了,要用统计学去估计。
⑷ 怎么利用FFT算法对音频噪声进行处理
fft是快速傅里叶变换,它是频谱分析的一种重要工具,例如,在处理过程中使用了快速傅立叶
变换fft,因此用平均周期图法计算功率谱密度函数估计是非常迅速的
⑸ 一维实序列的快速傅里叶变换(FFT)
通过前面的分析,我们认识到傅里叶变换本身是复数运算,地球物理获取的数据大多数是实数,对于实数的变换原则上可直接套用复序列的FFT算法,但那样是把实数序列当作虚部为零的复数对待,显然需要存储虚部的零并进行无功的运算,既浪费了一倍的计算内存,又降低了约一半的运算速度。
为了不浪费不可不设的虚部内存和必然出现的复数运算,可否将一个实序列分为两个子实序列,分别作为实部与虚部构成一个复数序列,然后用复序列的FFT算法求其频谱,对合成的复序列频谱进行分离和加工得到原实序列的频谱呢?答案是肯定的,实现这一过程思路就是实序列FFT算法的基本思想。
1.实序列的傅里叶变换性质
对于一个N个样本的实序列x(k),其频谱为X(j),用Xr(j)和Xi(j)表示X(j)的实部和虚部, 表示X(j)的共轭,则
证明:已知 则
地球物理数据处理基础
上式两端取共轭,并注意到x(k)是实序列,则
地球物理数据处理基础
这就是实序列的傅里叶变换具有复共轭性。
其同样具有周期性,即
地球物理数据处理基础
2.一维实序列的FFT算法
(1)同时计算两个实序列的FFT算法
已知两个实序列h(k),g(k)(k=0,1,…,N-1),例如重磁异常平面数据中的两条剖面,或地震勘探中的两道地震记录,可以人为地构成一个复序列:
地球物理数据处理基础
设h(k)的频谱为H(j)=Hr(j)+iHi(j)
g(k)的频谱为G(j)=Gr(j)+iGi(j)
y(k)的频谱为Y(j)=Yr(j)+i Yi(j)
利用上节的复序列FFT算法,求得Y(j),即Yr(j)和Yi(j)已知,来寻找Hr(j),Hi(j),Gr(j),Gi(j)与Yr(j),Yi(j)之间的关系。
对式(8-22)作傅里叶变换:
地球物理数据处理基础
由于H(j),G(j)本身是复序列,所以不能仅从上式分离出H(j)和G(j)。应用Y(j)的周期性,容易得到
Y(N-j)=H(-j)+iG(-j)
上式取共轭:
地球物理数据处理基础
由于h(k),g(k)为实序列,对上式右端应用复共轭定理,得
地球物理数据处理基础
对式(8-23)展开,得
地球物理数据处理基础
对式(8-24)展开,并应用共轭关系,得
地球物理数据处理基础
把式(8-25)和式(8-26)与Y(j)=Yr(j)+iYi(j)进行对比,有
地球物理数据处理基础
整理得
地球物理数据处理基础
因此,对于两个实序列,通过构造一个复序列,应用复序列的FFT算法和式(8-28)的分离加工,即可得到两个实序列的频谱。
(2)计算2 N个数据点的实序列FFT算法
设有2N点的实序列u(k)(k=0,1,…,2N-1),首先按k的偶、奇分成两个子实序列,并构成复序列,即
地球物理数据处理基础
通过调用复序列FFT算法,求得y(k)的频谱为Y(j)。另记h(k),g(k)的频谱为H(j)和G(j)。
利用前面式(8-23)和式(8-24),容易求得
地球物理数据处理基础
下面分析用H(j),G(j)形成u(k)频谱的问题。记u(k)(k=0,1,…,2 N-1)的频谱为V(j),分析V(j),H(j),G(j)之间的关系,根据定义
地球物理数据处理基础
利用式(8-31)和式(8-34)可换算出u(k)的前N个频谱V(j)(j=0,1,…,N-1),还要设法求u(k)的后N个频谱V(N+j)(j=0,1,…,N-1)。利用实序列其频谱的复共轭和周期性:
(1)H(N)=H(0),G(N)=G(0),WN1=-1,得
地球物理数据处理基础
(2)由于u(k)(k=0,1,…,2N-1)是实序列,同样利用实序列其频谱的复共轭和周期性,用已求出的前N个频谱V(j)表示出后面的N-1个频谱V(N+j):
地球物理数据处理基础
由于0<2N-j<N,所以可从V(j)(j=0,1,…,N-1)中选出V(2N-j)(j=N+1,N+2,…,2 N-1),并直接取其共轭 即可得到V(N+1)~V(2 N-1),从而完成整个实序列频谱的计算。
总结以上叙述,一维实序列u(k)(k=0,1,…,2N-1)的FFT计算编程步骤如下:
(1)按偶、奇拆分实序列u(k),并构造复序列:
地球物理数据处理基础
(2)调用复序列的FFT计算y(k)的频谱Y(j)(j=0,1,…,N-1);
(3)用下式计算形成h(k),g(k)的频谱H(j)和G(j);
地球物理数据处理基础
(4)用下式换算实序列u(k)的频谱V(j)(j=0,1,…,2 N-1):
地球物理数据处理基础
[例3]求实序列样本u(k)={1,2,1,1,3,2,1,2}(k=0,1,…,7)的频谱。
解:按偶、奇拆分实序列u(k),按式(8-37)构造复序列c(j)(j=0,1,2,3),即
c(0)=1+2i; c(1)=1+i; c(2)=3+2i; c(3)=1+2i。
(1)调用复序列FFT求c(j)(j=0,1,2,3)的频谱Z(k)(k=0,1,2,3),得
Z(0)=6+7i; Z(1)=-3; Z(2)=2+i; Z(3)=-1。
地球物理数据处理基础
(3)运用公式(8-38)计算H(j),G(j):
地球物理数据处理基础
(4)根据式(8-39)求出u(k)(k=0,1,…,7)的8个频谱V(j)(j=0,1,…,7):
地球物理数据处理基础
地球物理数据处理基础
由上例可见,完成全部2 N个实序列频谱的计算只需做N次FFT计算,相比直接用复序列的FFT算法节省了约一半的计算量。
⑹ 示波器怎么用fft实现频谱分析
具有此功能的示波器,点击MATH键,在屏上找到FFT功能并按此键即可。
⑺ matlab如何用fft
matlab自带的fft函数是快速傅里叶变换函数。主要用于降噪处理,通过使用傅里叶变换求噪声中隐藏的信号的频率分量。
该函数使用方法:
方法一:
Y= fft(X)用快速傅里叶变换 (FFT) 算法计算X的离散傅里叶变换(DFT)。
如果X是向量,则fft(X)返回该向量的傅里叶变换。
如果X是矩阵,则fft(X)将X的各列视为向量,并返回每列的傅里叶变换。
如果X是一个多维数组,则fft(X)将沿大小不等于 1 的第一个数组维度的值视为向量,并返回每个向量的傅里叶变换。
方法二:
如果X是向量且X的长度小于n,则为X补上尾零以达到长度n。
如果X是向量且X的长度大于n,则对X进行截断以达到长度n。
如果X是矩阵,则每列的处理与在向量情况下相同。
如果X为多维数组,则大小不等于 1 的第一个数组维度的处理与在向量情况下相同。
Y= fft(X,n)返回n点 DFT。如果未指定任何值,则Y的大小与X相同。
我们通过下例,来了解fft函数使用过程:
第一步、指定信号的参数,采样频率为 1 kHz,信号持续时间为 1.5 秒。
Fs=1000;%采样频率
T=1/Fs;%采样周期
L=1500;%信号长度
t=(0:L-1)*T;%时间向量
第二步、构造一个信号,其中包含幅值为 0.7 的 50 Hz 正弦量和幅值为 1 的 120 Hz 正弦量。
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
第三步、用均值为零、方差为 4 的白噪声扰乱该信号。
X = S + 2*randn(size(t));
第四步、在时域中绘制含噪信号。通过查看信号 X(t) 很难确定频率分量。
plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)'),ylabel('X(t)')
第五步、计算信号的傅里叶变换。
Y = fft(X);
第六步、计算双侧频谱 P2, 计算单侧频谱 P1。
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1)
第七步、定义频域 f 并绘制单侧幅值频谱 P1
f = Fs*(0:(L/2))/L;
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)'),ylabel('|P1(f)|')
运行结果。
⑻ 如何使用MATLAB中的fft函数来进行频谱分析
因为对于大多数实信号来说,正负频率是对称的,所以fft计算频谱后,实际是计算了在一个周期内的谱,而习惯上是只要在[0,Ts/2]观察就可以了,俗称单功率谱或者单边频谱,因此程序就取2了。你自己也可以全部取,注意x向量同样对应长度就可以了。
⑼ 如何用FFT得到谐波幅值频率和相位
用FFT得到谐波的频谱,里面含有频率,幅度和相位,同时可以通过这个三个而求得其他参数。
FFT是一种DFT的高效算法,称为快速傅立叶变换(fast Fourier transform),它根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。FFT对傅氏变换的理论并没有新的发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换,可以说是进了一大步。
⑽ 正弦序列FFT频谱分析程序问题!!
因为N个样点的信号经过fft以后变成N个样点的频谱,这个频谱是关于第N/2+1样点左右对称的,所以真正有用的频谱数据只有前面一半,后面一半是镜像。mxk11是对前N/2个样点取幅度谱,其实应该是取1:N1/2+1,你这里少取了一个点。具体为什么会镜像请看数字信号处理DFT章节。