基2演算法和直接演算法

① 請用MATLAB、C語言或者其他語言編程實現8點序列的基2-DIT－FFT演算法，並對結果進行分析驗證。

typedef struct

{

float real;

float img;

}COMPLEX;

void reverse(COMPLEX *X,COMPLEX *x,int N) /*輸入數組名稱和元素個數，實現倒序 */

{

int LH = N/2;

int j = LH;

int N1 = N-2;

int i;

for(i = 1; i <= N1; i++)

{

int k = LH;

X[i] = (i<j ? x[j] : x[i]);

/*

if(i < j)

{

comp T = a[j];

a[j] = a[i];

a[i] = T;

*/

}

while(j >= k)

{

j = j-k;

k = k/2;

}

j = j+k;

}

}

/*
函數功能：對指定長度的采樣數據進行快速傅立葉變換，並返回變換值, X = FFT(x)

入口參數：目標采樣數據，保存變換值的指針，采樣數據長度

出口參數：指向變換值的指針

注意事項：保存結果的空間是在函數外部申請，本函數不處理此內容

/*/

COMPLEX *FFT(COMPLEX *X,COMPLEX *x, unsigned N)

{

unsigned temp=N,L=0,M=0,LE=1,LE1=0,J=0,I=0,IP=0;

COMPLEX U={0,0},W={0,0},T={0,0};

while((temp>>=1)>0) /*獲得階碼M*/

{

M++;

}

reverse(X,x,N);

for(L=0;L<M;L++)

{

LE=1<<L; /*LE=2^L,分組間隔*/

LE1 = LE>>1; /*LE1=LE/2，對偶跨距，實際也是一個分組中對偶點的對數*/

U.real = 1.0;

U.img = 1.0;

W.real = cos(PI / LE1);

W.img = sin(PI / LE1);

J = 0;

while(J<LE1)

{

I = J;

while(I < N)

{

IP = I + LE1;

ComMul(X + IP,&U,&T);

ComSub(X + I,&T,X + IP);

ComAdd(X + I,&T,X + I);

I += LE;

}

ComMul(&U,&W,&U);

J++;

}

}

return X;

}

FFT中有幾個復數運算，自己實現，不想發給你，年輕人還是要自己動手做點東西。N為任意數，正常應該為2的冪次方。

② 誰知道DFT和FFT的發展歷史啊

DFT/FFT的發展歷史
離散傅里葉變換（Discrete Fourier Transform，DFT）是數字信號處理最重要的基石之一，也是對信號進行分析和處理時最常用的工具之一。在200多年前法國數學家、物理學家傅里葉提出後來以他名字命名的傅里葉級數之後，用DFT這個工具來分析信號就已經為人們所知。歷史上最偉大的數學家之一。
歐拉是第一個使用「函數」一詞來描述包含各種參數的表達式的人，例如：y = f(x)。他是把微積分應用於物理學的先驅者之一。 給出了一個用實變數函數表示傅立葉級數系數的方程； 用三角級數來描述離散聲音在彈性媒介中傳播，發現某些函數可以通過餘弦函數之和來表達。 但在很長時間內，這種分析方法並沒有引起更多的重視，最主要的原因在於這種方法運算量比較大。直到1965年，Cooley和Tukey在《計算機科學 》發表著名的《機器計算傅立葉級數的一種演算法》論文，FFT才開始大規模應用。
那個年代，有個肯尼迪總統科學咨詢委員會。其中有項研究主題是，對蘇聯核測試進行檢測，Tukey就是其中一員。美國/蘇聯核測試提案的批准，主要取決於不實地訪問核測試設施而做出檢測的方法的發展。其中一個想法是，分析離海岸的地震計情況，這種計算需要快速演算法來計算DFT。其它應用是國家安全，如用聲學探測遠距離的核潛艇。所以在軍事上，迫切需要一種快速的傅立葉變換演算法，這也促進了FFT的正式提出。
FFT的這種方法充分利用了DFT運算中的對稱性和周期性，從而將DFT運算量從N2減少到N*log2N。當N比較小時，FFT優勢並不明顯。但當N大於32開始，點數越大，FFT對運算量的改善越明顯。比如當N為1024時，FFT的運算效率比DFT提高了100倍。在庫利和圖基提出的FFT演算法中，其基本原理是先將一個N點時域序列的DFT分解為N個1點序列的DFT，然後將這樣計算出來的N個1點序列DFT的結果進行組合，得到最初的N點時域序列的DFT值。實際上，這種基本的思想很早就由德國偉大的數學家高斯提出過，在某種情況下，天文學計算（也是現在FFT應用的領域之一）與等距觀察的有限集中的行星軌道的內插值有關。由於當時計算都是靠手工，所以產生一種快速演算法的迫切需要。 而且，更少的計算量同時也代表著錯誤的機會更少，正確性更高。高斯發現，一個富氏級數有寬度N=N1*N2，可以分成幾個部分。計算N2子樣本DFT的N1長度和N1子樣本DFT的N2長度。只是由於當時尚欠東風——計算機還沒發明。在20世紀60年代，伴隨著計算機的發展和成熟，庫利和圖基的成果掀起了數字信號處理的革命，因而FFT發明者的桂冠才落在他們頭上。
之後，桑德（G.Sand）-圖基等快速演算法相繼出現，幾經改進，很快形成了一套高效運算方法，這就是現在的快速傅立葉變換（FFT）。這種演算法使DFT的運算效率提高1到2個數量級，為數字信號處理技術應用於各種信號的實時處理創造了良好的條件，大大推進了數學信號處理技術。1984年，法國的杜哈梅（P.Dohamel）和霍爾曼（H.Hollamann）提出的分裂基塊快速演算法，使運算效率進一步提高。
庫利和圖基的FFT演算法的最基本運算為蝶形運算，每個蝶形運算包括兩個輸入點，因而也稱為基-2演算法。在這之後，又有一些新的演算法，進一步提高了FFT的運算效率，比如基-4演算法，分裂基演算法等。這些新演算法對FFT運算效率的提高一般在50%以內，遠遠不如FFT對DFT運算的提高幅度。從這個意義上說，FFT演算法是里程碑式的。可以說，正是計算機技術的發展和FFT的出現，才使得數字信號處理迎來了一個嶄新的時代。除了運算效率的大幅度提高外，FFT還大大降低了DFT運算帶來的累計量化誤差，這點常為人們所忽略。

分給我吧 哈哈

③ 試畫出信號按時間抽取的基-2FFT演算法信號流圖，並寫出對應過程奇偶分流的公式。

基2演算法，序列的長度是為2的冪，序列的DFT為。序列可以由奇序列和偶序列組成，DFT分別為和。 從最後一級往前分解對應的蝶形結構，這些蝶形結構最左邊的輸入都是序列的DFT值，而分解直到最左邊的蝶形結構是兩點序列的DFT，此時最左邊的值是序列x[k]。

f1=50; %10Hz

f2=100; %100Hz

%抽樣頻率

Fs=1000; %100Hz

%抽樣點數N

L=10;

N=2^L;

%抽樣脈沖序列

n = 0:N-1;

t = n./Fs;

% f2 一個周期的采樣數

M = floor(Fs/f2);

%被采樣信號

x = cos(2*pi*f1.*t)+sin(2*pi*f2.*t);

%采樣序列

subplot(311);

stem(t(1:2*M),x(1:2*M));

hold off;

%傅里葉變換

%根據有限長序列的離散傅里葉變換公式計算DFT

n = 0:N-1;

k = 0:N-1;

F = x * exp(-j*2*pi/N).^(n'*k);

subplot(312);

plot(n,abs(F));

subplot(313);

plot(k,angle(F));

(3)基2演算法和直接演算法擴展閱讀：

庫利－圖基快速傅里葉變換演算法是將序列長為N的DFT分區為兩個長為N/2的子序列的DFT，因此這一應用只適用於序列長度為2的冪的DFT計算，即基2-FFT。實際上，如同高斯和庫利與圖基都指出的那樣，庫利－圖基演算法也可以用於序列長度N為任意因數分解形式的DFT，即混合基FFT，而且還可以應用於其他諸如分裂基FFT等變種。

盡管庫利－圖基演算法的基本思路是採用遞歸的方法進行計算，大多數傳統的演算法實現都將顯示的遞歸演算法改寫為非遞歸的形式。另外，因為庫利－圖基演算法是將DFT分解為較小長度的多個DFT，因此它可以同任一種其他的DFT演算法聯合使用。

④ 基-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綆楁硶鐨勫熀鏈鎬濊礬鏄浠涔堬紵

鍩2綆楁硶錛屽簭鍒楃殑闀垮害鏄涓2鐨勫籙錛屽簭鍒楃殑DFT涓恆傚簭鍒楀彲浠ョ敱濂囧簭鍒楀拰鍋跺簭鍒楃粍鎴愶紝DFT鍒嗗埆涓哄拰銆 浠庢渶鍚庝竴綰у線鍓嶅垎瑙ｅ瑰簲鐨勮澏褰㈢粨鏋勶紝榪欎簺銛跺艦緇撴瀯鏈宸﹁竟鐨勮緭鍏ラ兘鏄搴忓垪鐨凞FT鍊礆紝鑰屽垎瑙ｇ洿鍒版渶宸﹁竟鐨勮澏褰㈢粨鏋勬槸涓ょ偣搴忓垪鐨凞FT錛屾ゆ椂鏈宸﹁竟鐨勫兼槸搴忓垪x[k]銆

鍩4鏃墮棿鎶藉彇FFT璁＄畻錛氬皢搴忓垪鍒嗕負4涓鐭搴忓垪錛屽垎鍒涓簒[4k]銆亁[4k+1]銆亁[4k+2]銆亁[4k+3]錛屾瘡涓綰ф湁N/4涓銛跺艦榪愮畻錛岀涓綰ф瘡涓銛跺艦榪愮畻涓嶉渶瑕佷箻閫夋嫨鍥犲瓙錛屾墍浠ユ病鏈夊嶆暟涔樻硶銆備箣鍚庣殑綰ф暟錛屽洜涓猴紝涓嶉渶瑕佷箻縐錛屾墍浠ユ瘡涓銛跺艦榪愮畻閮介渶瑕3嬈′箻鏃嬭漿鍥犲瓙錛屾晠闇瑕3嬈″嶆暟涔樻硶銆

鎵╁睍璧勬枡錛

娉ㄦ剰浜嬮」錛

1銆佹渶灝忚垂鐢ㄦ渶澶ф祦鍩轟簬鏈澶ф祦闂棰橈紝閭ｄ箞灝辮佹敞鎰忔祦鐨勪笁涓閲嶈佺殑鎬ц川錛屽叾嬈¤佺悊瑙ｆ祦閲忔槸涓縐嶉熺巼,鑰屼笉鏄鎬婚噺銆

2銆佹渶澶ф祦闂棰樼殑涓縐嶈В娉曞氨鏄澧炲箍璺錛屾渶閲嶈佺殑灝辨槸瑕佺悊瑙ｉ嫻佹搷浣滐紝鏄涓烘挙閿鍏堝墠閫夊彇鐨勪笉鍚堥傜殑寮с

3銆佹瘡涓緇撶偣瑕佸緩絝4涓寮э紝鍒嗕負涓ゅ癸紝姣忓歸兘鏄涓涓瀹歸噺涓鴻竟瀹歸噺鐨勬ｈ垂鐢╟ 鍜 鍙嶅悜瀹歸噺涓0鐨勪粯璐圭敤錛屼袱瀵圭殑鍖哄埆鏄鏂瑰悜鎯沖弽錛屾瘡瀵圭殑鎿嶄綔鐢ㄥ紓鎴朸銆

鍙傝冭祫鏂欐潵婧愶細鐧懼害鐧劇-FFT鍘熺悊