軟體fft演算法

❶ IDFT的運算怎麼用FFT實現，在軟體中怎樣進行設置

IDFT就是Inverse Discrete Fourier Transform 離散傅里葉逆變換。

FFT就是Fast Fourier Transform 快速傅里葉變換。

兩者的應用都是將時域中難以處理的信號轉換成易於復處理的頻域信號，分析完成後進行傅里葉反變換即得到原始的粗瞎鄭時域岩頌信號。

兩者的異同是：我們知道在數學上用級數來無限逼進某個函數，以便簡化計算過程而又不致使誤差過大，這樣工程上才能應用，否則一些數學模型是無法實現快速求解的。

IDFT：對於有限長的序列我們可以使用離散傅立葉變換，神指IDFT是對制序列傅立葉變換的等距采樣。

FFT：並不是與IDFT不相同的另一種變換（即原理是一樣的），而是為了減少IDFT運算次數的一種快速演算法。它是對IDFT變換式進行一次次的分解，使其成為若干小點數IDFT的組合，從而減小運算量。常用的FFT是以2為基數，它的運算效率高，程序比較簡單，使用也十分地方便。

❷ 怎樣用C語言實現FFT演算法啊

1、二維FFT相當於對行和列分別進行一維FFT運算。具體的實現辦法如下：
先對各行逐一進行一維FFT，然後再對變換後的新矩陣的各列逐一進行一維FFT。相應的偽代碼如下所示：
for (int i=0; i<M; i++)
FFT_1D(ROW[i]，N);
for (int j=0; j<N; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩陣的第i行。注意這只是一個簡單的記法，並不能完全照抄。還需要通過一些語句來生成各行的數據。同理，COL[i]是對矩陣的第i列的一種簡單表示方法。
所以，關鍵是一維FFT演算法的實現。

2、常式：

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#defineN1000
/*定義復數類型*/
typedefstruct{
doublereal;
doubleimg;
}complex;
complexx[N],*W;/*輸入序列,變換核*/
intsize_x=0;/*輸入序列的大小，在本程序中僅限2的次冪*/
doublePI;/*圓周率*/
voidfft();/*快速傅里葉變換*/
voidinitW();/*初始化變換核*/
voidchange();/*變址*/
voidadd(complex,complex,complex*);/*復數加法*/
voidmul(complex,complex,complex*);/*復數乘法*/
voidsub(complex,complex,complex*);/*復數減法*/
voidoutput();
intmain(){
inti;/*輸出結果*/
system("cls");
PI=atan(1)*4;
printf("Pleaseinputthesizeofx:
");
scanf("%d",&size_x);
printf("Pleaseinputthedatainx[N]:
");
for(i=0;i<size_x;i++)
scanf("%lf%lf",&x[i].real,&x[i].img);
initW();
fft();
output();
return0;
}
/*快速傅里葉變換*/
voidfft(){
inti=0,j=0,k=0,l=0;
complexup,down,proct;
change();
for(i=0;i<log(size_x)/log(2);i++){/*一級蝶形運算*/
l=1<<i;
for(j=0;j<size_x;j+=2*l){/*一組蝶形運算*/
for(k=0;k<l;k++){/*一個蝶形運算*/
mul(x[j+k+l],W[size_x*k/2/l],&proct);
add(x[j+k],proct,&up);
sub(x[j+k],proct,&down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*初始化變換核*/
voidinitW(){
inti;
W=(complex*)malloc(sizeof(complex)*size_x);
for(i=0;i<size_x;i++){
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
}
}
/*變址計算，將x(n)碼位倒置*/
voidchange(){
complextemp;
unsignedshorti=0,j=0,k=0;
doublet;
for(i=0;i<size_x;i++){
k=i;j=0;
t=(log(size_x)/log(2));
while((t--)>0){
j=j<<1;
j|=(k&1);
k=k>>1;
}
if(j>i){
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
/*輸出傅里葉變換的結果*/
voidoutput(){
inti;
printf("Theresultareasfollows
");
for(i=0;i<size_x;i++){
printf("%.4f",x[i].real);
if(x[i].img>=0.0001)printf("+%.4fj
",x[i].img);
elseif(fabs(x[i].img)<0.0001)printf("
");
elseprintf("%.4fj
",x[i].img);
}
}
voidadd(complexa,complexb,complex*c){
c->real=a.real+b.real;
c->img=a.img+b.img;
}
voidmul(complexa,complexb,complex*c){
c->real=a.real*b.real-a.img*b.img;
c->img=a.real*b.img+a.img*b.real;
}
voidsub(complexa,complexb,complex*c){
c->real=a.real-b.real;
c->img=a.img-b.img;
}

❸ 如何應用matlab進行fft分析

FFT是離散傅立葉變換的快速演算法，可以將一個信號變換
到頻域。有些信號在時域上是很難看出什麼特徵的，但是如
果變換到頻域之後，就很容易看出特徵了。這就是很多信號
分析採用FFT變換的原因。另外，FFT可以將一個信號的頻譜
提取出來，這在頻譜分析方面也是經常用的。
雖然很多人都知道FFT是什麼，可以用來做什麼，怎麼去
做，但是卻不知道FFT之後的結果是什意思、如何決定要使用
多少點來做FFT。

現在圈圈就根據實際經驗來說說FFT結果的具體物理意義。
一個模擬信號，經過ADC采樣之後，就變成了數字信號。采樣
定理告訴我們，采樣頻率要大於信號頻率的兩倍，這些我就
不在此羅嗦了。

采樣得到的數字信號，就可以做FFT變換了。N個采樣點，
經過FFT之後，就可以得到N個點的FFT結果。為了方便進行FFT
運算，通常N取2的整數次方。

假設采樣頻率為Fs，信號頻率F，采樣點數為N。那麼FFT
之後結果就是一個為N點的復數。每一個點就對應著一個頻率
點。這個點的模值，就是該頻率值下的幅度特性。具體跟原始
信號的幅度有什麼關系呢？假設原始信號的峰值為A，那麼FFT
的結果的每個點（除了第一個點直流分量之外）的模值就是A
的N/2倍。而第一個點就是直流分量，它的模值就是直流分量
的N倍。而每個點的相位呢，就是在該頻率下的信號的相位。
第一個點表示直流分量（即0Hz），而最後一個點N的再下一個
點（實際上這個點是不存在的，這里是假設的第N+1個點，也
可以看做是將第一個點分做兩半分，另一半移到最後）則表示
采樣頻率Fs，這中間被N-1個點平均分成N等份，每個點的頻率
依次增加。例如某點n所表示的頻率為：Fn=(n-1)*Fs/N。
由上面的公式可以看出，Fn所能分辨到頻率為為Fs/N，如果
采樣頻率Fs為1024Hz，采樣點數為1024點，則可以分辨到1Hz。
1024Hz的采樣率采樣1024點，剛好是1秒，也就是說，采樣1秒
時間的信號並做FFT，則結果可以分析到1Hz，如果采樣2秒時
間的信號並做FFT，則結果可以分析到0.5Hz。如果要提高頻率
分辨力，則必須增加采樣點數，也即采樣時間。頻率解析度和
采樣時間是倒數關系。
假設FFT之後某點n用復數a+bi表示，那麼這個復數的模就是
An=根號a*a+b*b，相位就是Pn=atan2(b,a)。根據以上的結果，
就可以計算出n點（n≠1，且n<=N/2）對應的信號的表達式為：
An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。
對於n=1點的信號，是直流分量，幅度即為A1/N。
由於FFT結果的對稱性，通常我們只使用前半部分的結果，
即小於采樣頻率一半的結果。

好了，說了半天，看著公式也暈，下面圈圈以一個實際的
信號來做說明。

假設我們有一個信號，它含有2V的直流分量，頻率為50Hz、
相位為-30度、幅度為3V的交流信號，以及一個頻率為75Hz、
相位為90度、幅度為1.5V的交流信號。用數學表達式就是如下：

S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)

式中cos參數為弧度，所以-30度和90度要分別換算成弧度。
我們以256Hz的采樣率對這個信號進行采樣，總共采樣256點。
按照我們上面的分析，Fn=(n-1)*Fs/N，我們可以知道，每兩個
點之間的間距就是1Hz，第n個點的頻率就是n-1。我們的信號
有3個頻率：0Hz、50Hz、75Hz，應該分別在第1個點、第51個點、
第76個點上出現峰值，其它各點應該接近0。實際情況如何呢？
我們來看看FFT的結果的模值如圖所示。

圖1 FFT結果
從圖中我們可以看到，在第1點、第51點、和第76點附近有
比較大的值。我們分別將這三個點附近的數據拿上來細看：
1點： 512+0i
2點： -2.6195E-14 - 1.4162E-13i
3點： -2.8586E-14 - 1.1898E-13i

50點：-6.2076E-13 - 2.1713E-12i
51點：332.55 - 192i
52點：-1.6707E-12 - 1.5241E-12i

75點：-2.2199E-13 -1.0076E-12i
76點：3.4315E-12 + 192i
77點：-3.0263E-14 +7.5609E-13i

很明顯，1點、51點、76點的值都比較大，它附近的點值
都很小，可以認為是0，即在那些頻率點上的信號幅度為0。
接著，我們來計算各點的幅度值。分別計算這三個點的模值，
結果如下：
1點： 512
51點：384
76點：192
按照公式，可以計算出直流分量為：512/N=512/256=2；
50Hz信號的幅度為：384/(N/2)=384/(256/2)=3；75Hz信號的
幅度為192/(N/2)=192/(256/2)=1.5。可見，從頻譜分析出來
的幅度是正確的。
然後再來計算相位信息。直流信號沒有相位可言，不用管
它。先計算50Hz信號的相位，atan2(-192, 332.55)=-0.5236,
結果是弧度，換算為角度就是180*(-0.5236)/pi=-30.0001。再
計算75Hz信號的相位，atan2(192, 3.4315E-12)=1.5708弧度，
換算成角度就是180*1.5708/pi=90.0002。可見，相位也是對的。
根據FFT結果以及上面的分析計算，我們就可以寫出信號的表達
式了，它就是我們開始提供的信號。

總結：假設采樣頻率為Fs，采樣點數為N，做FFT之後，某
一點n（n從1開始）表示的頻率為：Fn=(n-1)*Fs/N；該點的模值
除以N/2就是對應該頻率下的信號的幅度（對於直流信號是除以
N）；該點的相位即是對應該頻率下的信號的相位。相位的計算
可用函數atan2(b,a)計算。atan2(b,a)是求坐標為(a,b)點的角
度值，范圍從-pi到pi。要精確到xHz，則需要采樣長度為1/x秒
的信號，並做FFT。要提高頻率解析度，就需要增加采樣點數，
這在一些實際的應用中是不現實的，需要在較短的時間內完成
分析。解決這個問題的方法有頻率細分法，比較簡單的方法是
采樣比較短時間的信號，然後在後面補充一定數量的0，使其長度
達到需要的點數，再做FFT，這在一定程度上能夠提高頻率分辨力。
具體的頻率細分法可參考相關文獻。

[附錄：本測試數據使用的matlab程序]
close all; %先關閉所有圖片
Adc=2; %直流分量幅度
A1=3; %頻率F1信號的幅度
A2=1.5; %頻率F2信號的幅度
F1=50; %信號1頻率(Hz)
F2=75; %信號2頻率(Hz)
Fs=256; %采樣頻率(Hz)
P1=-30; %信號1相位(度)
P2=90; %信號相位(度)
N=256; %采樣點數
t=[0:1/Fs:N/Fs]; %采樣時刻

%信號
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
%顯示原始信號
plot(S);
title('原始信號');

figure;
Y = fft(S,N); %做FFT變換
Ayy = (abs(Y)); %取模
plot(Ayy(1:N)); %顯示原始的FFT模值結果
title('FFT 模值');

figure;
Ayy=Ayy/(N/2); %換算成實際的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; %換算成實際的頻率值
plot(F(1:N/2),Ayy(1:N/2)); %顯示換算後的FFT模值結果
title('幅度-頻率曲線圖');

figure;
Pyy=[1:N/2];
for i="1:N/2"
Pyy(i)=phase(Y(i)); %計算相位
Pyy(i)=Pyy(i)*180/pi; %換算為角度
end;
plot(F(1:N/2),Pyy(1:N/2)); %顯示相點陣圖
title('相位-頻率曲線圖');

看完這個你就明白諧波分析了。

❹ 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的公式是什麼和演算法是怎樣實現

二維FFT相當於對行和列分別進行一維FFT運算。具體的實現辦法如下：
先對各行逐一進行一維FFT，然後再對變換後的新矩陣的各列逐一進行一維FFT。相應的偽代碼如下所示：
for (int i=0; i<M; i++)
FFT_1D(ROW[i]，N);
for (int j=0; j<N; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩陣的第i行。注意這只是一個簡單的記法，並不能完全照抄。還需要通過一些語句來生成各行的數據。同理，COL[i]是對矩陣的第i列的一種簡單表示方法。
所以，關鍵是一維FFT演算法的實現。下面討論一維FFT的演算法原理。

【1D-FFT的演算法實現】
設序列h(n)長度為N，將其按下標的奇偶性分成兩組，即he和ho序列，它們的長度都是N/2。這樣，可以將h(n)的FFT計算公式改寫如下 ：

(A)
由於

所以，(A)式可以改寫成下面的形式：

按照FFT的定義，上面的式子實際上是：

其中，k的取值范圍是 0~N-1。
我們注意到He(k)和Ho(k)是N/2點的DFT，其周期是N/2。因此，H(k)DFT的前N/2點和後N/2點都可以用He(k)和Ho(k)來表示

❻ Python科學計算——復雜信號FFT

FFT (Fast Fourier Transform, 快速傅里葉變換) 是離散傅里葉變換的快速演算法，也是數字信號處理技術中經常會提到的一個概念。用快速傅里葉變換能將時域的數字信號轉換為頻域信號，轉換為頻域信號後我們可以很方便地分析出信號的頻率成分。

當我們把雙頻信號FFT示例中的 fft_size 的值改為 2**12 時，這時，基頻為 16Hz，不能被 1kHz整除，所以 1kHz 處發生了頻譜泄露，而它能被 4kHz 整除，所以 4kHz 可以很好地被采樣。

由於波形的前後不是連續的，出現波形跳變，而跳變處有著非常廣泛的頻譜，因此FFT的結果中出現了頻譜泄漏。

為了減小FFT所截取的數據段前後的跳變，可以對數據先乘以一個窗函數，使得其前後數據能平滑過渡。常用的hanning窗函數的定義如下：

50Hz 正弦波與hann窗函數乘積之後的重復波形如下：

我們對頻譜泄漏示例中的1kHz 和 4kHz 信號進行了 hann 窗函數處理，可以看出能量更加集中在 1kHz 和 4kHz，在一定程度上抑制了頻譜泄漏。

以 1kHz 三角波為例，我們知道三角波信號中含有豐富的頻率信息，它的傅里葉級數展開為：

當數字信號的頻率隨時間變化時，我們稱之為掃頻信號。以頻率隨時間線性變化的掃頻信號為例，其數學形式如下：

其頻率隨時間線性變化，當我們在 [0,1] 的時間窗口對其進行采樣時，其頻率范圍為 0~5kHz。當時間是連續時，掃頻信號的頻率也是連續的。但是在實際的處理中，是離散的點采樣，因此時間是不連續的，這就使掃頻信號的快速傅里葉變換問題退化為多點頻信號快速傅里葉變換問題。其快速傅里葉變換得到的頻譜圖如下所示：

以 50Hz 正弦信號相位調制到 1kHz 的信號為例，其信號形式如下：

它的時域波形，頻率響應和相位響應如下圖所示：

以掃頻信號為例，當我們要探究FFT中的能量守恆時，我們要回歸到信號最初的形式：

❼ 實序列的FFT演算法

在以上討論FFT演算法中，均假定序列x(l)為復的，但實際問題中的序列大多為實的。當然，我們可以把實序列處理成虛部為零的復序列。因此，就要引進許多零參加運算。這樣一來，在機器運算時間和存儲單元方面都將造成很大的浪費。在本段中，我們介紹對實序列x(l)應用FFT演算法的一個有效方法。

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

設有N=4的兩個實序列x₁(l)與x₂(l)。為了求得它們的譜X₁(m)與X₂(m)，我們用此二實序列構造成如下復序列

物探數字信號分析與處理技術

利用上一段的方法，可以求得復序列x(l)的譜X(m)。根據(7-3-1)得到

物探數字信號分析與處理技術

上式中的m用N-m代替，則得

物探數字信號分析與處理技術

將上式兩端取共軛，根據對稱性有

物探數字信號分析與處理技術

根據DFT的復共軛性質，對於實序列x₁(l)與x₂(l)，有

物探數字信號分析與處理技術

於是從(7-3-4)得到

物探數字信號分析與處理技術

聯立求解(7-3-2)和(7-3-6)便得到

物探數字信號分析與處理技術

例如設有兩個N=4點的實序列，

物探數字信號分析與處理技術

我們用它們構造一個N=4點的復序列

物探數字信號分析與處理技術

利用FFT演算法求X(m)，m=0，1，2，3(圖7-3-1)，

圖7-3-1 N=4點的FFT演算法流程圖

於是得到

物探數字信號分析與處理技術

因此從式(7-3-7)得到

物探數字信號分析與處理技術

物探數字信號分析與處理技術

2.實序列的FFT演算法

設有N點的實序列x(l)，l=0，1，2，…，N-1。按照點的奇偶編號，將它們分成N/2個點的兩個子序列

物探數字信號分析與處理技術

設x₁(l)的譜與x₂(l)的譜分別為X₁(m)與X₂(m)

物探數字信號分析與處理技術

其中

於是可以將實序列x(l)的譜X(m)，用兩個子序列x₁(l)，x₂(l)的譜X₁(m)，X₂(m)來表示

物探數字信號分析與處理技術

其中

物探數字信號分析與處理技術

注意，x₁(l)，x₂(l)與X₁(m)，X₂(m)均以N/2為周期，

利用x₁(l)、x₂(l)構成如下復序列

物探數字信號分析與處理技術

利用FFT演算法可以求得復序列 的譜 。根據(7-3-7)就求得兩個實子序列的譜X₁(m)與X₂(m)

物探數字信號分析與處理技術

有了X₁(m)，X₂(m)，根據(7-3-10)就可求得X(m)。以上就是用FFT演算法求實序列x(l)的譜X(m)的方法。必須指出，用公式(7-3-10)求X(m)時，第一，兩個實子序列的譜X₁(m)，X₂(m)及復序列x珓(l)的譜珘X(m)均是以N/2為周期的周期序列;第二，由於x

(l)是實序列，根據DFT的復共軛性質有X(m)=X^*(N-m)，m=0，1，…，N/2，故只需求得前(N/2)+1個點的X(m)，就得到全部N個點的X(m)了

例如，有N=8點的實序列，

物探數字信號分析與處理技術

首先，按點的奇偶編號分成兩個實子序列，

物探數字信號分析與處理技術

其次用它們構造如下復序列，

物探數字信號分析與處理技術

用FFT演算法求此復序列的譜 (圖7-3-2)

圖7-3-2 N=4點的FFT演算法流程圖

於是得到:

根據周期性，有

物探數字信號分析與處理技術

根據(7-3-12)式，

物探數字信號分析與處理技術

根據周期性，有

物探數字信號分析與處理技術

故最終由(7-3-10)得到

物探數字信號分析與處理技術

❽ 16點DFT的FFT演算法

FFT（快速傅里葉變換）是DFT的一種特殊情況，就是當運算點的個數是2的整數次冪的時候進行的運算（不夠用0補齊）。

FFT計算原理及流程圖：

原理：FFT的計算要求點數必須為2的整數次冪，如果點數不夠用0補齊。例如計算{2，3，5，8，4}的16點FFT，需要補11個0後進行計算。FFT計算運用蝶形運算，在蝶形運算中變化規律由W(N, p)推導，其中N為FFT計算點數，J為下角標的值。

L = 1時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0;

L = 2時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1;

L = 3時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1, 2, 3;

所以，W(N, p) = W(2^L, J)，其中J = 0, 1, ..., 2^(L-1)-1

又因為2^L = 2^M*2^(L-M) = N*2^(L-M)，這里N為2的整數次冪，即N=2^M，

W(N, p) = W(2^L, J) = W(N*2^(L-M), J) = W(N, J*2^(M-L))

所以，p = J*2^(M-L)，此處J = 0, 1, ..., 2^(L-1)-1，當J遍歷結束但計算點數不夠N時，J=J+2^L，後繼續遍歷，直到計算點數為N時不再循環。

流程圖：

/*======================================================================
*方法名：fft
*方法功能：計算數組的FFT，運用蝶形運算
*
*變數名稱：
*yVector-原始數據
*length-原始數據長度
*N-FFT計算點數
*fftYreal-FFT後的實部
*fftYImg-FFT後的虛部
*
*返回值：是否成功的標志，若成功返回true，否則返回false
*=====================================================================*/

+(BOOL)fft:(floatfloat*)yVectorandOriginalLength:(NSInteger)lengthandFFTCount:(NSInteger)NandFFTReal:(floatfloat*)fftYRealandFFTYImg:(floatfloat*)fftYImg
{
//確保計算時時2的整數冪點數計算
NSIntegerN1=[selfnextNumOfPow2:N];

//定義FFT運算是否成功的標志
BOOLisFFTOK=false;

//判斷計算點數是否為2的整數次冪
if(N!=N1)
{
//不是2的整數次冪，直接計算DFT
isFFTOK=[selfdft:yVectorandOriginalLength:lengthandFFTCount:NandFFTReal:fftYRealandFFTYImg:fftYImg];

//返回成功標志
returnisFFTOK;
}


//如果計算點數位2的整數次冪，用FFT計算，如下
//定義變數
floatyVectorN[N1];//N點運算的原始數據
NSIntegerpowOfN=log2(N1);//N=2^powOfN，用於標記最大運算級數(公式中表示為：M)
NSIntegerlevel=1;//運算級數（第幾次運算），最大為powOfN，初值為第一級運算(公式中表示為：L)
NSIntegerlineNum;//行號，倒序排列後的蝶形運算行號(公式中表示為：k)
floatinverseOrderY[N1];//yVector倒序x
NSIntegerdistanceLine=1;//行間距，第level級運算每個蝶形的兩個節點距離為distanceLine=2^(L-1)(公式中表示為：B)
NSIntegerp;//旋轉因子的階數，旋轉因子表示為W(N,p)，p=J*2^(M-L)
NSIntegerJ;//旋轉因子的階數，旋轉因子表示為W(2^L,J)，J=0,1,2,...,2^(L-1)-1=distanceLine-1
floatrealTemp,imgTemp,twiddleReal,twiddleImg,twiddleTheta,twiddleTemp=PI_x_2/N1;
NSIntegerN_4=N1/4;

//判斷點數是否夠FFT運算點數
if(length<=N1)
{
//如果N至少為length，先把yVector全部賦值
for(NSIntegeri=0;i<length;i++)
{
yVectorN[i]=yVector[i];
}

if(length<N1)
{
//如果N>length後面補零
for(NSIntegeri=length;i<N1;i++)
{
yVectorN[i]=0.0;
}
}
}
else
{
//如果N<length截取相應長度的數據進行運算
for(NSIntegeri=0;i<N1;i++)
{
yVectorN[i]=yVector[i];
}
}

//調用倒序方法
[selfinverseOrder:yVectorNandN:N1andInverseOrderVector:inverseOrderY];

//初始值
for(NSIntegeri=0;i<N1;i++)
{
fftYReal[i]=inverseOrderY[i];
fftYImg[i]=0.0;
}

//三層循環
//第三層（最里）：完成相同旋轉因子的蝶形運算
//第二層（中間）：完成旋轉因子的變化，步進為2^level
//第一層（最外）：完成M次迭代過程，即計算出x(k)=A0(k),A1(k),...,Am(k)=X(k)

//第一層循環
while(level<=powOfN)
{
distanceLine=powf(2,level-1);//初始條件distanceLine=2^(level-1)
J=0;
NSIntegerpow2_Level=distanceLine*2;//2^level
NSIntegerpow2_NSubL=N1/pow2_Level;//2^(M-L)

//第二層循環
while(J<distanceLine)
{
p=J*pow2_NSubL;
lineNum=J;
NSIntegerstepCount=0;//J運算的步進計數

//求旋轉因子
if(p==0)
{
twiddleReal=1.0;
twiddleImg=0.0;
}
elseif(p==N_4)
{
twiddleReal=0.0;
twiddleImg=-1.0;
}
else
{
//計算尤拉公式中的θ
twiddleTheta=twiddleTemp*p;

//計算復數的實部與虛部
twiddleReal=cos(twiddleTheta);
twiddleImg=-11*sin(twiddleTheta);
}

//第三層循環
while(lineNum<N1)
{
//計算下角標
NSIntegerfootNum=lineNum+distanceLine;

/*---------------------------------------
*用復數運算計算每級中各行的蝶形運算結果
*X(k)=X(k)+X(k+B)*W(N,p)
*X(k+B)=X(k)-X(k+B)*W(N,p)
*---------------------------------------*/
realTemp=fftYReal[footNum]*twiddleReal-fftYImg[footNum]*twiddleImg;
imgTemp=fftYReal[footNum]*twiddleImg+fftYImg[footNum]*twiddleReal;

//將計算後的實部和虛部分別存放在返回數組中
fftYReal[footNum]=fftYReal[lineNum]-realTemp;
fftYImg[footNum]=fftYImg[lineNum]-imgTemp;
fftYReal[lineNum]=fftYReal[lineNum]+realTemp;
fftYImg[lineNum]=fftYImg[lineNum]+imgTemp;

stepCount+=pow2_Level;

//行號改變
lineNum=J+stepCount;
}

//旋轉因子的階數變換，達到旋轉因子改變的效果
J++;
}

//運算級數加一
level++;
}

isFFTOK=true;
returnisFFTOK;
}

❾ 基-2fft演算法的軟體實現 matlab代碼

% 基於Matlab的時間抽取基2FFT演算法
function y=myditfft(x)
%本程序對輸入序列實現DIT-FFT基2演算法，點數取大於等於長度的2的冪次
%------------------------------------
% Leo's fft program(改編網上的一個程序)
%------------------------------------
m=log2(2^nextpow2(length(x))); %求的x長度對應的2的最低冪次m
N=2^m;
if length(x)<N
x=[x,zeros(1,N-length(x))]; %若長度不是2的冪，補0到2的整數冪
end
x;
%--------------------------------------------------------------------------
%對輸入序列進行倒序
%如果輸入序列的自然順序號I用二進制數（例如n2n1n0）表示
%則其倒位序J對應的二進制數就是（n0n1n2），這樣，在原來自然順序時應該放x(I)的
%單元，現在倒位序後應放x(J)。
%--------------------------------------------------------------------------
%以下程序相當於以下程序:
%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1：2^m數列的倒序
%y=x(nxd); %將倒序排列作為初始值
%--------------------------------------------------------------------------
NV2=N/2;
NM1=N-1;
I=0;
J=0;
while I<NM1
if I<J
T=x(J+1);
x(J+1)=x(I+1);
x(I+1)=T;
end
K=NV2;

while K<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
x;
%--------------------------------------------------------------------------
%以下程序解釋:
%第一級從x(0)開始,跨接一階蝶形,再取每條對稱
%第二級從x(0)開始,跨接兩階蝶形,再取每條對稱
%第m級從x(0)開始,跨接2^(m-1)階蝶形,再取每條對稱....
%--------------------------------------------------------------------------
for mm=1:m %將DFT做m次基2分解，從左到右，對每次分解作DFT運算
Nmr=2^mm;
u=1; %旋轉因子u初始化
WN=exp(-j*2*pi/Nmr); %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
for n=1:Nmr/2 %本次跨越間隔內的各次碟形運算
for k=n:Nmr:N %本次碟形運算的跨越間隔為Nmr=2^mm
kp=k+Nmr/2; %確定碟形運算的對應單元下標(對稱性)
t=x(kp)*u; %碟形運算的乘積項
x(kp)=x(k)-t; %碟形運算的加法項
x(k)=x(k)+t;
end
u=u*WN; %修改旋轉因子，多乘一個基本DFT因子WN
end
end
y=x; %輸出

❿ FFT演算法分幾種

FFT演算法分析FFT演算法的基本原理是把長序列的DFT逐次分解為較短序列的DFT。按照抽取方式的不同可分為DIT-FFT（按時間抽取）和DIF-FFT（按頻率抽取）演算法。按照蝶形運算的構成不同可分為基2、基4、基8以及任意因子（2n,n為大於1的整數），基2、基4演算法較為常用。 網上有幫助文檔: http://www.5doc.com/doc/123035(右上角有點擊下載)