fft頻譜演算法

⑴ 二維實序列的快速傅里葉變換（FFT）

在地球物理數據處理中，經常遇到處理二維實數據的情況。例如在地震勘探中，對面波勘探數據作頻散分析解釋時，要將時間－空間域的信息轉換為頻率－波數域頻譜；在重磁異常的濾波或轉換中，要將空間域的異常f（x，y）轉換為波數域F（ω，υ）等。這些分析都需要進行二維的傅里葉變換（FFT）。

根據傅里葉變換的定義，對於連續二維函數f（x，y），其傅里葉變換對為

地球物理數據處理基礎

對於離散的二維序列f_jk（j＝0，1，…，M－1；k＝0，1，…，N－1），其傅里葉變換為

地球物理數據處理基礎

1.二維復序列的FFT演算法

對於M條測線，每條測線N個測點，構成復序列y_jk（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）將y_jk輸入A₁（k），即A₁（k）＝y_jk（k＝0，1，…，N－1）；

3）計算W^r，存入W（r），即

4）q＝1，2，…，p（p＝log₂N），若q為偶數執行6），否則執行5）；

5）k＝0，1，2，…，（2^p－q－1）和n＝0，1，2，…，（2^q－1－1）循環，作

A₂（k2^q＋n）＝A₁（k2^q－1＋n）＋A₁（k2^q－1＋n＋2^p－1）

A₂（k2^q＋n＋2^q－1）＝［A₁（k2^q－1＋n）－A₁（k2^q－1＋n＋2^p－1）］·W（k2^q－1）

至k，n循環結束；

6）k＝0，1，2，…，（2^p－q－1）和n＝0，1，2，…，（2^q－1－1）循環，作

A₁（k2^q＋n）＝A₂（k2^q－1＋n）＋A₂（k2^q－1＋n＋2^p－1）

A₁（k2^q＋n＋2^q－1）＝［A₂（k2^q－1＋n）－A₂（k2^q－1＋n＋2^p－1）］·W（k2^q－1）

至k，n循環結束；

7）q循環結束，若p為偶數，將A₁（n）輸入到Y_jn，否則將A₂（n）輸入到Y_jn（n＝0，1，…，N－1）；

8）j循環結束，得到Y_jn（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）將Y_jn輸入A₁（j），即A₁（j）＝Y_jn（j＝0，1，…，M－1）；

3）計算W^r存入W（r），即

4）q＝1，2，…，p（p＝log₂M），若q為偶數執行6），否則執行5）；

5）j＝0，1，2，…，（2^p－q－1）和m＝0，1，2，…，（2^q－1－1）循環，作

A₂（j2^q＋m）＝A₁（j2^q－1＋m）＋A₁（j2^q－1＋m＋2^p－1）

A₂（j2^q＋m＋2^q－1）＝［A₁（j2^q－1＋m）－A₁（j2^q－1＋m＋2^p－1）］·W（j2^q－1）

至j，m循環結束；

6）j＝0，1，2，…，（2^p－q－1）和m＝0，1，2，…，（2^q－1－1）循環，作

A₁（j2^q＋m）＝A₂（j2^q－1＋m）＋A₂（j2^q－1＋m＋2^p－1）

A₁（j2^q＋m＋2^q－1）＝［A₂（j2^q－1＋m）－A₂（j2^q－1＋m＋2^p－1）］·W（j2^q－1）

至j，m循環結束；

7）q循環結束，若p為偶數，將A₁（m）輸入到Y_mn，否則將A₂（m）輸入到Y_mn（m＝0，1，…，M－1）；

8）n循環結束，得到二維復序列的傅氏變換Y_mn（m＝0，1，…，M－1；n＝0，1，…，N－1），

所求得的Y_mn是復數值，可以寫為

Y_mn＝R_mn＋iI_mn （m＝0，1，…，M－1；n＝0，1，…，N－1）

其中，R_mn，I_mn的值也是已知的。

2.二維實序列的FFT演算法

對於二維的實序列，我們把其看作是虛部為零的復序列，套用上述的二維復序列FFT方法來求其頻譜演算法上也是可行的，但勢必會增加大量的無功運算。因此，有必要研究二維實序列FFT的實用演算法，同一維實序列FFT的實現思路一樣，同樣把二維實序列按一定的規律構造成二維復序列，調用二維復序列FFT，然後通過分離和加工得到原實序列的頻譜。

例如采樣區域有2 M條測線，每條測線有N個點，並且M，N都是2的整數冪，需要計算實樣本序列x_jk（j＝0，1，2，…，2 M－1；k＝0，1，2，…，N－1）的傅氏變換：

地球物理數據處理基礎

類似於一維實序列FFT的思想，直接建立下面的二維實序列FFT演算法：

（1）將一個二維實序列按偶、奇線號分為兩個二維子實序列，分別作為實部和虛部組合為一個二維復序列。即令

地球物理數據處理基礎

（2）調用二維復FFT過程，求出y_jk的二維傅氏變換Y_mn的復數值：

地球物理數據處理基礎

式中：R_mn，I_mn是Y_mn的實部和虛部。

（3）利用R_mn，I_mn換算X_mn的值。

前兩步容易實現，下面分析第（3）步的實現。

記h_jk，g_jk的傅氏變換為H_mn，G_mn。根據傅里葉變換的定義，我們導出X_mn與H_mn，G_mn的關系式：

地球物理數據處理基礎

式中，H_mn，G_mn為復數，我們用上標r和i表示其實部和虛部，將上式右端實部、虛部分離

地球物理數據處理基礎

其中：

地球物理數據處理基礎

下面的任務是將H_mn，G_mn各分量與通過二維復FFT求出的R_mn，I_mn值聯系起來。為此先給出奇、偶分解性質和類似於一維情況的三個二維傅氏變換性質：

（1）奇偶分解性

任何一個正負對稱區間定義的函數，均可唯一地分解為如下偶（even）、奇（odd）函數之和：

地球物理數據處理基礎

（2）周期性

地球物理數據處理基礎

（3）復共軛性

地球物理數據處理基礎

現在我們來建立R_mn，I_mn與H_mn，G_mn的關系。對Y_mn作奇偶分解：

地球物理數據處理基礎

根據線性性質

地球物理數據處理基礎

對照式（8－54）和式（8－55），得

地球物理數據處理基礎

由於h_jk，g_jk是實函數，根據復共軛性質，上面兩式對應的奇偶函數相等。即

地球物理數據處理基礎

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

地球物理數據處理基礎

將式（8－57）代入式（8－50），得

地球物理數據處理基礎

再利用H_mn，G_mn周期性及復共軛性，可以得到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）即可得到原始序列x_jk（j＝0，1，…，2M－1；n＝0，1，…，N－1）在m＝0，1，…，M－1；n＝0，1，…，N－1區間的傅氏變換X_mn。

具體二維實序列的FFT演算法如下：

（1）令h_jk＝x_2j，k，g_jk＝x_2j＋1，k，形成

y_jk＝h_jk＋ig_jk （j＝0，1，…，2 M－1；n＝0，1，…，N－1）

（2）調用二維復序列FFT過程，即從兩個方向先後調用一維復FFT演算法式（8－43）和式（8－44），求得y_jk的二維傅氏變換Y_mn的復數值：

Y_mn＝R_mn＋iI_mn （m＝0，1，…，M－1；n＝0，1，…，N－1）

（3）用下列公式由R_mn，I_mn的值換算X_mn的值：

地球物理數據處理基礎

地球物理數據處理基礎

⑵ 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），用X_r（j）和X_i（j）表示X（j）的實部和虛部， 表示X（j）的共軛，則

證明：已知 則

地球物理數據處理基礎

上式兩端取共軛，並注意到x（k）是實序列，則

地球物理數據處理基礎

這就是實序列的傅里葉變換具有復共軛性。

其同樣具有周期性，即

地球物理數據處理基礎

2.一維實序列的FFT演算法

（1）同時計算兩個實序列的FFT演算法

已知兩個實序列h（k），g（k）（k＝0，1，…，N－1），例如重磁異常平面數據中的兩條剖面，或地震勘探中的兩道地震記錄，可以人為地構成一個復序列：

地球物理數據處理基礎

設h（k）的頻譜為H（j）＝H_r（j）＋iH_i（j）

g（k）的頻譜為G（j）＝G_r（j）＋iG_i（j）

y（k）的頻譜為Y（j）＝Y_r（j）＋i Y_i（j）

利用上節的復序列FFT演算法，求得Y（j），即Y_r（j）和Y_i（j）已知，來尋找H_r（j），H_i（j），G_r（j），G_i（j）與Y_r（j），Y_i（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）＝Y_r（j）＋iY_i（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），W^N₁＝－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 的第一個數組維度的值視為向量，並返回每個向量的傅里葉變換。

• 方法二：

Y= fft(X,n)返回n點 DFT。如果未指定任何值，則Y的大小與X相同。

• 如果X是向量且X的長度小於n，則為X補上尾零以達到長度n。

• 如果X是向量且X的長度大於n，則對X進行截斷以達到長度n。

• 如果X是矩陣，則每列的處理與在向量情況下相同。

• 如果X為多維數組，則大小不等於 1 的第一個數組維度的處理與在向量情況下相同。

我們通過下例，來了解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章節。