查表實現的FFT演算法運行時間

1. FFT的演算法

FFT是一種DFT的高效演算法，稱為快速傅立葉變換（fast Fourier transform），它根據離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的演算法進行改進獲得的。FFT演算法可分為按時間抽取演算法和按頻率抽取演算法，先簡要介紹FFT的基本原理。從DFT運算開始，說明FFT的基本原理。
DFT的運算為：

式中
由這種方法計算DFT對於 的每個K值，需要進行4N次實數相乘和（4N-2）次相加，對於N個k值，共需4N*4N次實數相乘和（4N-2)(4N-2）次實數相加。改進DFT演算法，減小它的運算量，利用DFT中 的周期性和對稱性，使整個DFT的計算變成一系列迭代運算，可大幅度提高運算過程和運算量，這就是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的優越性。

2. TMS320F2812的FFT32點，64點，128點，256點的計算時間多長阿

以TI的DSP2812測試，系統時鍾150M，程序調入RAM中運行。
分別測試512點FFT，256點FFT和128點FFT四個大步驟用時。
512點 用時(us)
1步倒序(RFFT32_brev) 42.4
2步計算(fft.cal) 223.8
3步實部虛部分離(fft.split) 13.1
4步計算幅值(fft.mag) 13.8
總用時 293.1
256點 用時(us)
1步倒序(RFFT32_brev) 21.5
2步計算(fft.cal) 94.9
3步實部虛部分離(fft.split) 13.1
4步計算幅值(fft.mag) 13.8
總用時 143.3
128點 用時(us)
1步倒序(RFFT32_brev) 10.8
2步計算(fft.cal) 39
3步實部虛部分離(fft.split) 13.1
4步計算幅值(fft.mag) 13.8
總用時 76.7

3. 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;
}

4. 我做「基於FFT演算法與實現」和「FIR濾波器的設計與實現」的實驗。。

1.1 實驗目的

1．了解數字信號處理系統的一般構成；

2．掌握奈奎斯特抽樣定理。

1.2 實驗儀器

1．YBLD智能綜合信號源測試儀 1台

2．雙蹤示波器 1台

3．MCOM－TG305數字信號處理與現代通信技術實驗箱 1台

4．PC機（裝有MATLAB、MCOM－TG305配套實驗軟體） 1台

1.3 實驗原理

一個典型的DSP系統除了數字信號處理部分外，還包括A/D和D/A兩部分。這是因為自然界的信號，如聲音、圖像等大多是模擬信號，因此需要將其數字化後進行數字信號處理，模擬信號的數字化即稱為A/D轉換。數字信號處理後的數據可能需還原為模擬信號，這就需要進行D/A轉換。一個僅包括A/D和D/A兩部分的簡化數字信號處理系統功能如圖1所示。

A/D轉換包括三個緊密相關的過程，即抽樣、量化和編碼。A/D轉換中需解決的以下幾個重要問題：抽樣後輸出信號中還有沒有原始信號的信息？如果有能不能把它取出來？抽樣頻率應該如何選擇？

奈奎斯特抽樣定理（即低通信號的均勻抽樣定理）告訴我們，一個頻帶限制在0至fx以內的低通信號x(t)，如果以fs≥2fx的抽樣速率進行均勻抽樣，則x(t)可以由抽樣後的信號xs(t)完全地確定，即xs(t)包含有x(t)的成分，可以通過適當的低通濾波器不失真地恢復出x(t)。最小抽樣速率fs=2fx稱為奈奎斯特速率。

低通

解碼

編碼

量化

抽樣

輸入信號 樣點輸出 濾波輸出

A/D（模數轉換） D/A（數模轉換）

圖1 低通采樣定理演示

為方便實現，實驗中更換了一種表現形式，即抽樣頻率固定（10KHz），通過改變輸入模擬信號的頻率來展示低通抽樣定理。我們可以通過研究抽樣頻率和模擬信號最高頻率分量的頻率之間的關系，來驗證低通抽樣定理。

1.4 實驗內容

1．軟體模擬實驗：編寫並調試MATLAB程序，分析有關參數，記錄有關波形。

2．硬體實驗：輸入不同頻率的正弦信號，觀察采樣時鍾波形、輸入信號波形、樣點輸出波形和濾波輸出波形。

1.5 MATLAB參考程序和模擬內容

%*******************************************************************%

%f—餘弦信號的頻率

% M—基2 FFT冪次數 N=2^M為采樣點數，這樣取值是為了便於作基2的FFT分析

%2. 采樣頻率Fs

%*******************************************************************%

function samples(f,Fs,M)

N=2^M; % fft點數=取樣總點數

Ts=1/Fs; % 取樣時間間隔

T=N*Ts; % 取樣總時間=取樣總點數*取樣時間間隔

n=0:N-1;

t=n*Ts;

Xn=cos(2*f*pi*t);

subplot(2,1,1);

stem(t,Xn);

axis([0 T 1.1*min(Xn) 1.1*max(Xn)]);

xlabel('t -->');

ylabel('Xn');

Xk=abs(fft(Xn,N));

subplot(2,1,2);

stem(n,Xk);

axis([0 N 1.1*min(Xk) 1.1*max(Xk)]);

xlabel('frequency -->');

ylabel('!Xk!');

%*******************************************************************%

假如有一個1Hz的餘弦信號y=cos(2*π*t)，對其用4Hz的采樣頻率進行采樣，共采樣32點，只需執行samples(1,4,5)，即可得到模擬結果。

軟體模擬實驗內容如下表所示：

模擬參數
f
Fs
Wo(計算)
Xn（圖形）
Xk（圖形）

（1，4，5）

另外記錄圖形，並標圖號

（1，8，5）

（2，8，6）

自 選

1.6 硬體實驗步驟

本實驗箱采樣頻率fs固定為10KHz，低通濾波器的截止頻率約為4.5KHz。

1、用低頻信號源產生正弦信號，正弦信號源頻率f自定，並將其接至2TP2（模擬輸入）端，將示波器通道一探頭接至2TP6（采樣時鍾）端觀察采樣時鍾波形，示波器通道二探頭接至2TP2觀察並記錄輸入信號波形。

2、將示波器通道二探頭接至2TP3觀察並記錄樣點輸出波形。

3、將示波器通道二探頭接至2TP4觀察並記錄濾波輸出波形。

4、根據采樣定理，分f=fs /8、f=fs/4、f=fs/2等3種情況更改正弦信號頻率，重復步驟2至步驟3。

5、用低頻信號源產生方波信號，重復步驟1至步驟4。

1.7 思考題

1、 討論在模擬實驗中所計算的數字域頻率Wo和Xk的圖形中非零譜線位置之間的對應關系。

2、 討論在模擬實驗中自選參數的意義。

3、將在2TP2端加方波信號後的恢復波形，與相同頻率的正弦信號的恢復波形相比，能夠得出哪些結論？

2 FFT頻譜分析實驗

2.1 實驗目的

1．通過實驗加深對快速傅立葉變換（FFT）基本原理的理解。

2．了解FFT點數與頻譜解析度的關系，以及兩種加長序列FFT與原序列FFT的關系。

2.2 實驗儀器

1．YBLD智能綜合信號源測試儀 1台

2．雙蹤示波器 1台

3．MCOM－TG305數字信號處理與現代通信技術實驗箱 1台

4．PC機（裝有MATLAB、MCOM－TG305配套實驗軟體） 1台

2.3 實驗原理

離散傅里葉變換（DFT）和卷積是信號處理中兩個最基本也是最常用的運算，它們涉及到信號與系統的分析與綜合這一廣泛的信號處理領域。實際上卷積與DFT之間有著互通的聯系：卷積可化為DFT來實現，其它的許多演算法，如相關、濾波和譜估計等都可化為DFT來實現，DFT也可化為卷積來實現。

對N點序列x(n)，其DFT變換對定義為：

在DFT運算中包含大量的重復運算。FFT演算法利用了蝶形因子WN的周期性和對稱性，從而加快了運算的速度。FFT演算法將長序列的DFT分解為短序列的DFT。N點的DFT先分解為2個N/2點的DFT，每個N/2點的DFT又分解為2個N/4點的DFT。按照此規律，最小變換的點數即所謂的「基數（radix）。」因此，基數為2的FFT演算法的最小變換（或稱蝶形）是2點DFT。一般地，對N點FFT，對應於N個輸入樣值，有N個頻域樣值與之對應。一般而言，FFT演算法可以分為時間抽取（DIT）FFT和頻率抽取（DIF）兩大類。

在實際計算中，可以採用在原來序列後面補0的加長方法來提高FFT的解析度；可以採用在原來序列後面重復的加長方法來增加FFT的幅度。

2.4 實驗內容

1．軟體模擬實驗：分別觀察並記錄正弦序列、方波序列及改變FFT的點數後的頻譜；分別觀察並記錄正弦序列、方波序列及2種加長序列等信號的頻譜。

2．硬體實驗：分別觀察並記錄正弦信號、方波信號及改變FFT的點數後的頻譜。

2.5 MATLAB參考程序和模擬內容

%*******************************************************************%

function[x]=ffts(mode,M)

Nfft=2^M;

x=zeros(1,Nfft); %定義一個長度為Nfft的一維全0數組

if mode= =1 for n=0:Nfft-1 x(n+1)=sin(2*pi*n/Nfft); end

end %定義一個長度為Nfft的單周期正弦序列

if mode= =2 for n=0:Nfft-1 x(n+1)=sin(4*pi*n/Nfft); end

end %定義一個長度為Nfft的雙周期正弦序列

if mode= =3 for n=0:Nfft/2-1 x(n+1)=sin(4*pi*n/Nfft); end

end %定義一個長度為Nfft/2的正弦序列，後面一半為0序列。

if mode= =4 for n=0:Nfft-1 x(n+1)=square(2*pi*n/Nfft); end

end

if mode= =5 for n=0:Nfft-1 x(n+1)=square(2*pi*n/Nfft); end

end

if mode= =6 for n=0:Nfft/2-1 x(n+1)=square(4*pi*n/Nfft); end

end

n=0:Nfft-1;

subplot(2,1,1);

stem(n,x);

axis([0 Nfft-1 1.1*min(x) 1.1*max(x)]);

xlabel('Points-->');

ylabel('x(n)');

X=abs(fft(x,Nfft));

subplot(2,1,2);

stem(n,X);

axis([0 Nfft-1 1.1*min(X) 1.1*max(X)]);

xlabel('frequency-->');

ylabel('!X(k)!');

%*******************************************************************%

假設需觀察方波信號的頻譜，對一個周期的方波信號作32點的FFT，則只需在MATLAB的命令窗口下鍵入：[x]=ffts(21,5) ，程序進行模擬，並且輸出FFT的結果。

關於軟體模擬實驗內容，建議在完成大量模擬例子的基礎上，選擇能夠體現實驗要求的4個以上的例子進行記錄。例如要觀察後面補0的加長方法來提高FFT的解析度的現象，可以模擬ffts(4,5)和ffts(6,6)兩個例子。

2.6 硬體實驗步驟

1．將低頻信號源輸出加到實驗箱模擬通道1輸入端，將示波器探頭接至模擬通道1輸出端。

2．在保證實驗箱正確加電且串口電纜連接正常的情況下，運行數字信號處理與DSP應用實驗開發軟體，在「數字信號處理實驗」菜單下選擇「FFT頻譜分析」子菜單，出現顯示FFT頻譜分析功能提示信息的窗口。

3．用低頻信號產生器產生一個1KHz的正弦信號。

4．選擇FFT頻譜分析與顯示的點數為64點，開始進行FFT運算。此後，計算機將周期性地取回DSP運算後的FFT數據並繪圖顯示

5．改信號源頻率，觀察並記錄頻譜圖的變化。

6．選擇FFT的點數為128點，觀察並記錄頻譜圖的變化。

7．更改正弦信號的頻率，重復步驟4 ~步驟6。

8．用低頻信號產生器產生一個1KHz的方波信號，重復步驟4 ~步驟7。注意：應根據實驗箱采樣頻率fs為10KHz和方波信號的頻帶寬度選擇方波信號的頻率。

本硬體實驗要進行兩種信號，每個信號兩種頻率，每個信號兩種點數等共8次具體實驗內容，性質能夠體現實驗要求的4個以上的例子進行記錄。

2.7 思考題

1．對同一個信號，不同點數FFT觀察到的頻譜圖有何區別？

2．序列加長後FFT與原序列FFT的關系是什麼，試推導其中一種關系。

3．用傅立葉級數理論，試說明正弦信號頻譜和方波信號頻譜之間的關系。

3 IIR濾波器設計實驗

3.1 實驗目的

1．通過實驗加深對IIR濾波器基本原理的理解。

2．學習編寫IIR濾波器的MATLAB模擬程序。

3.2 實驗儀器

1．YBLD智能綜合信號源測試儀 1台

2．雙蹤示波器 1台

3．MCOM－TG305數字信號處理與現代通信技術實驗箱 1台

4．PC機（裝有MATLAB、MCOM－TG305配套實驗軟體） 1台

3.3 實驗原理

IIR濾波器有以下幾個特點：

1．IIR數字濾波器的系統函數可以寫成封閉函數的形式。

2．IIR數字濾波器採用遞歸型結構，即結構上帶有反饋環路。IIR濾波器運算結構通常由延時、乘以系數和相加等基本運算組成，可以組合成直接型、正准型、級聯型、並聯型四種結構形式，都具有反饋迴路。由於運算中的舍入處理，使誤差不斷累積，有時會產生微弱的寄生振盪。

3．IIR數字濾波器在設計上可以藉助成熟的模擬濾波器的成果，如巴特沃斯、契比雪夫和橢圓濾波器等，有現成的設計數據或圖表可查，其設計工作量比較小，對計算工具的要求不高。在設計一個IIR數字濾波器時，我們根據指標先寫出模擬濾波器的公式，然後通過一定的變換，將模擬濾波器的公式轉換成數字濾波器的公式。

4．IIR數字濾波器的相位特性不好控制，對相位要求較高時，需加相位校準網路。

在MATLAB下設計IIR濾波器可使用Butterworth函數設計出巴特沃斯濾波器，使用Cheby1函數設計出契比雪夫I型濾波器，使用Cheby2設計出契比雪夫II型濾波器，使用ellipord函數設計出橢圓濾波器。下面主要介紹前兩個函數的使用。

與FIR濾波器的設計不同，IIR濾波器設計時的階數不是由設計者指定，而是根據設計者輸入的各個濾波器參數（截止頻率、通帶濾紋、阻帶衰減等），由軟體設計出滿足這些參數的最低濾波器階數。在MATLAB下設計不同類型IIR濾波器均有與之對應的函數用於階數的選擇。

一、巴特沃斯IIR濾波器的設計

在MATLAB下，設計巴特沃斯IIR濾波器可使用butter函數。

Butter函數可設計低通、高通、帶通和帶阻的數字和模擬IIR濾波器，其特性為使通帶內的幅度響應最大限度地平坦，但同時損失截止頻率處的下降斜度。在期望通帶平滑的情況下，可使用butter函數。

butter函數的用法為：

[b,a]=butter(n,Wn,/ftype/)

其中n代表濾波器階數，Wn代表濾波器的截止頻率，這兩個參數可使用buttord函數來確定。buttord函數可在給定濾波器性能的情況下，求出巴特沃斯濾波器的最小階數n，同時給出對應的截止頻率Wn。buttord函數的用法為：

[n,Wn]= buttord(Wp,Ws,Rp,Rs)

其中Wp和Ws分別是通帶和阻帶的拐角頻率（截止頻率），其取值范圍為0至1之間。當其值為1時代表采樣頻率的一半。Rp和Rs分別是通帶和阻帶區的波紋系數。

不同類型（高通、低通、帶通和帶阻）濾波器對應的Wp和Ws值遵循以下規則：

1．高通濾波器：Wp和Ws為一元矢量且Wp>Ws；

2．低通濾波器：Wp和Ws為一元矢量且Wp<Ws；

3．帶通濾波器：Wp和Ws為二元矢量且Wp<Ws，如Wp=[0.2,0.7],Ws=[0.1,0.8];

4．帶阻濾波器：Wp和Ws為二元矢量且Wp>Ws，如Wp=[0.1,0.8],Ws=[0.2,0.7]。

二、契比雪夫I型IIR濾波器的設計

在期望通帶下降斜率大的場合，應使用橢圓濾波器或契比雪夫濾波器。在MATLAB下可使用cheby1函數設計出契比雪夫I型IIR濾波器。

cheby1函數可設計低通、高通、帶通和帶阻契比雪夫I型濾IIR波器，其通帶內為等波紋，阻帶內為單調。契比雪夫I型的下降斜度比II型大，但其代價是通帶內波紋較大。

cheby1函數的用法為：

[b,a]=cheby1(n,Rp,Wn,/ftype/)

在使用cheby1函數設計IIR濾波器之前，可使用cheblord函數求出濾波器階數n和截止頻率Wn。cheblord函數可在給定濾波器性能的情況下，選擇契比雪夫I型濾波器的最小階和截止頻率Wn。

cheblord函數的用法為：

[n,Wn]=cheblord(Wp,Ws,Rp,Rs)

其中Wp和Ws分別是通帶和阻帶的拐角頻率（截止頻率），其取值范圍為0至1之間。當其值為1時代表采樣頻率的一半。Rp和Rs分別是通帶和阻帶區的波紋系數。

3.4 實驗內容

1．軟體模擬實驗：編寫並調試MATLAB程序，選擇不同形式，不同類型的4種濾波器進行模擬，記錄幅頻和相頻特性，對比巴特沃斯濾波器和契比雪夫濾波器。

2．硬體實驗：設計IIR濾波器，在計算機上觀察沖激響應、幅頻特性和相頻特性，然後下載到實驗箱。用示波器觀察輸入輸出波形，測試濾波器的幅頻響應特性。

3.5 MATLAB參考程序和模擬內容

%*******************************************************************%

%mode: 1--巴特沃斯低通；2--巴特沃斯高通；3--巴特沃斯帶通；4--巴特沃斯帶阻

% 5--契比雪夫低通；6--契比雪夫高通；7--契比雪夫帶通；8--契比雪夫帶阻

%fp1,fp2： 通帶截止頻率，當高通或低通時只有fp1有效

%fs1, fs2： 阻帶截止頻率，當高通或低通時只有fs1有效

%rp: 通帶波紋系數

%as: 阻帶衰減系數

%sample: 采樣率

%h: 返回設計好的濾波器系數

%*******************************************************************%

function[b,a]=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample)

wp1=2*fp1/sample;wp2=2*fp2/sample;

ws1=2*fs1/sample;ws2=2*fs2/sample;

%得到巴特沃斯濾波器的最小階數N和3bd頻率wn

if mode<3[N,wn]=buttord(wp1,ws1,rp,as);

elseif mode<5[N,wn]=buttord([wp1 wp2],[ws1 ws2],rp,as);

%得到契比雪夫濾波器的最小階數N和3bd頻率wn

elseif mode<7[N,wn]=cheb1ord(wp1,ws1,rp,as);

else[N,wn]=cheblord([wp1 wp2],[ws1 ws2],rp,as);

end

%得到濾波器系數的分子b和分母a

if mode= =1[b,a]=butter(N,wn);end

if mode= =2[b,a]=butter(N,wn,/high/);end

if mode= =3[b,a]=butter(N,wn);end

if mode= =4[b,a]=butter(N,wn,/stop/);end

if mode= =5[b,a]=cheby1(N,rp,wn);end

if mode= =6[b,a]=cheby1(N,rp,wn,/high/);end

if mode= =7[b,a]=cheby1(N,rp,wn);end

if mode= =8[b,a]=cheby1(N,rp,wn,/stop/);end

set(gcf,/menubar/,menubar);

freq_response=freqz(b,a);

magnitude=20*log10(abs(freq_response));

m=0:511;

f=m*sample/(2*511);

subplot(3,1,1);plot(f,magnitude);grid; %幅頻特性

axis([0 sample/2 1.1*min(magnitude) 1.1*max(magnitude)]);

ylabel('Magnitude');xlabel('Frequency-->');

phase=angle(freq_response);

subplot(3,1,2);plot(f,phase);grid; %相頻特性

axis([0 sample/2 1.1*min(phase) 1.1*max(phase)]);

ylabel('Phase');xlabel('Frequency-->');

h=impz(b,a,32); %32點的單位函數響應

t=1:32;

subplot(3,1,3);stem(t,h);grid;

axis([0 32 1.2*min(h) 1.1*max(h)]);

ylabel('h(n)');xlabel('n-->');

%*******************************************************************%

假設需設計一個巴特沃斯低通IIR濾波器，通帶截止頻率為2KHz，阻帶截止頻率為3KHz，通帶波紋系數為1，阻帶衰減系數為20，采樣頻率為10KHz，則只需在MATLAB的命令窗口下鍵入：

[b,a]=iirfilt(1,2000,3000,2400,2600,1,20,10000)

程序進行模擬，並且按照如下順序輸出數字濾波器系統函數

的系數

b= b0 b1 ……bn

a= a0 a1 ……an

關於軟體模擬實驗內容，建議在完成大量模擬例子的基礎上，選擇能夠體現實驗要求的4個例子進行記錄，系統函數只要記錄系統的階數。

3.6 硬體實驗步驟

1．根據實驗箱采樣頻率fs為10KHz的條件，用低頻信號發生器產生一個頻率合適的低頻正弦信號，將其加到實驗箱模擬通道1輸入端，將示波器通道1探頭接至模擬通道1輸入端，通道2探頭接至模擬通道2輸出端。

2．在保證實驗箱正確加電且串口電纜連接正常的情況下，運行數字信號處理與DSP應用實驗開發軟體，在「數字信號處理實驗」菜單下選擇「IIR濾波器」子菜單，出現提示信息。

3．輸入濾波器類型、濾波器截止頻率等參數後，分別點擊「幅頻特性」和「相頻特性」按鈕，在窗口右側觀察IIR濾波器的幅頻特性和相頻特性。此時提示信息將消失，如需查看提示信息，可點擊「設計說明」按鈕。

4．點擊「下載實現」按鈕，IIR濾波器開始工作，此時窗口右側將顯示IIR濾波器的幅頻特性。

5．根據輸入濾波器類型，更改低頻信號源的頻率，觀察示波器上輸入輸出波形幅度的變化情況，測量IIR濾波器的幅頻響應特性，看其是否與設計的幅頻特性一致。

6．更改濾波器類型、濾波器截止頻率等參數（共4種），重復步驟3至步驟5。所選擇的例子參數最好和MATLAB模擬程序的例子一樣。

7．用低頻信號產生器產生一個500Hz的方波信號，分別設計3種濾波器，完成如下表要求的功能，並且記錄參數和波形。

功 能
濾波器類型
參 數
輸出波形

fp1
fp2
fs1
fs2

通過3次及以下次數的諧波

另外記錄圖形，並標圖號

濾除5次及以下次數的諧波

通過3次到5次的諧波

3.7 思考題

1．在實驗箱采樣頻率fs固定為10KHz的條件下，要觀察方波信號頻帶寬度內的各個諧波分量，方波信號的頻率最高不能超過多少，為什麼？

2．硬體實驗內容7中輸出信號各個諧波分量，與原來方波信號同樣諧波分量相比，有沒有發生失真？主要發生了什麼類型的失真？為什麼？

4 窗函數法FIR濾波器設計實驗

4.1 實驗目的

1．通過實驗加深對FIR濾波器基本原理的理解。

2．學習使用窗函數法設計FIR濾波器，了解窗函數的形式和長度對濾波器性能的影響。

4.2 實驗儀器

1．YBLD智能綜合信號源測試儀 1台

2．雙蹤示波器 1台

3．MCOM－TG305數字信號處理與現代通信技術實驗箱 1台

4．PC機（裝有MATLAB、MCOM－TG305配套實驗軟體） 1台

4.3 實驗原理

數字濾波器的設計是數字信號處理中的一個重要內容。數字濾波器設計包括FIR（有限單位脈沖響應）濾波器與IIR（無限單位脈沖響應）濾波器兩種。

與IIR濾波器相比，FIR濾波器在保證幅度特性滿足技術要求的同時，很容易做到嚴格的線性相位特性。設FIR濾波器單位脈沖響應h(n)長度為N，其系統函數H（z）為：

H（z）是z－1的N－1次多項式，它在z平面上有N－1個零點，原點z=0是N－1階重極點，因此H（z）是永遠穩定的。穩定和線性相位特性是FIR濾波器突出的優點。

FIR濾波器的設計任務是選擇有限長度的h(n)。使傳輸函數H( )滿足技術要求。FIR濾波器的設計方法有多種，如窗函數法、頻率采樣法及其它各種優化設計方法，本實驗介紹窗函數法的FIR濾波器設計。

窗函數法是使用矩形窗、三角窗、巴特利特窗、漢明窗、漢寧窗和布萊克曼窗等設計出標准響應的高通、低通、帶通和帶阻FIR濾波器。

一、firl函數的使用

在MATLAB下設計標准響應FIR濾波器可使用firl函數。firl函數以經典方法實現加窗線性相位FIR濾波器設計，它可以設計出標準的低通、帶通、高通和帶阻濾波器。firl函數的用法為：

b=firl(n,Wn,/ftype/,Window)

各個參數的含義如下：

b—濾波器系數。對於一個n階的FIR濾波器，其n+1個濾波器系數可表示為：b(z)=b(1)+b(2)z－1+…+b(n+1)z－n。

n—濾波器階數。

Wn—截止頻率，0≤Wn≤1，Wn=1對應於采樣頻率的一半。當設計帶通和帶阻濾波器時，Wn=[W1 W2],W1≤ω≤W2。

ftype—當指定ftype時，可設計高通和帶阻濾波器。Ftype=high時，設計高通FIR濾波器；ftype=stop時設計帶阻FIR濾波器。低通和帶通FIR濾波器無需輸入ftype參數。

Window—窗函數。窗函數的長度應等於FIR濾波器系數個數，即階數n+1。

二、窗函數的使用

在MATLAB下，這些窗函數分別為：

1．矩形窗：w=boxcar(n)，產生一個n點的矩形窗函數。

2．三角窗：w=triang(n)，產生一個n點的三角窗函數。

當n為奇數時，三角窗系數為w(k)=
當n為偶數時，三角窗系數為w(k)=
3．巴特利特窗：w=Bartlett(n)，產生一個n點的巴特利特窗函數。

巴特利特窗系數為w(k)=
巴特利特窗與三角窗非常相似。巴特利特窗在取樣點1和n上總以零結束，而三角窗在這些點上並不為零。實際上，當n為奇數時bartlett(n)的中心n－2個點等效於triang(n－2)。

4．漢明窗：w=hamming(n)，產生一個n點的漢明窗函數。

漢明窗系數為w(k+1)=0.54－0.46cos（ ） k=0,…,n－1

5．漢寧窗：w=hanning(n)，產生一個n點的漢寧窗函數。

漢寧窗系數為w(k)=0.5[1－cos( )] k=1,…,n

6．布萊克曼窗：w=Blackman(n),產生一個n點的布萊克曼窗函數。

布萊克曼窗系數為w(k)=0.42－0.5cos(2π )+0.8cos（4π ）] k=1,…,n

與等長度的漢明窗和漢寧窗相比，布萊克曼窗的主瓣稍寬，旁瓣稍低。

7．凱澤窗：w=Kaiser(n,beta)，產生一個n點的凱澤窗數，其中beta為影響窗函數旁瓣的β參數，其最小的旁瓣抑制α與β的關系為:

0.1102(α－0.87) α>50

β= 0.5842(α－21)0.4+0.07886(α－21) 21≤α≤50

0 α<21

增加β可使主瓣變寬，旁瓣的幅度降低。

8．契比雪夫窗：w=chebwin(n,r)產生一個n點的契比雪夫窗函數。其傅里葉變換後的旁瓣波紋低於主瓣r個db數。

4.4 實驗內容

1．軟體模擬實驗：編寫並調試MATLAB程序，觀察不同窗，不同類型濾波器不同點數等共4種FIR濾波器的h(n)，並記錄幅頻特性和相頻特性。

2．硬體實驗：用窗函數法設計標准響應的FIR濾波器，在計算機上觀察窗函數幅頻特性、幅頻特性和相頻特性，然後下載到實驗箱。用示波器觀察輸入輸出波形，測試濾波器的幅頻響應特性。

4.5 MATLAB參考程序和模擬內容

%*******************************************************************%

%mode: 模式（1--高通；2--低通；3--帶通；4--帶阻）

%n: 階數，加窗的點數為階數加1

%fp: 高通和低通時指示截止頻率，帶通和帶阻時指示下限頻率

%fs: 帶通和帶阻時指示上限頻率

%window:加窗（1--矩形窗；2--三角窗；3--巴特利特窗；4--漢明窗；

% 5--漢寧窗；6--布萊克曼窗；7--凱澤窗；8--契比雪夫窗）

%r: 代表加chebyshev窗的r值和加kaiser窗時的beta值

%sample: 采樣率

%h: 返回設計好的FIR濾波器系數

%*******************************************************************%

%mode: 模式（1--高通；2--低通；3--帶通；4--帶阻）

%n: 階數，加窗的點數為階數加1

%fp: 高通和低通時指示截止頻率，帶通和帶阻時指示下限頻率

%fs:

5. 求用C++實現的FFT演算法

很早以前的，如果管用別忘了給我加分呀
/*
This computes an in-place complex-to-complex FFT
x and y are the real and imaginary arrays of 2^m points.
dir = 1 gives forward transform
dir = -1 gives reverse transform
*/
short FFT(short int dir,long m,double *x,double *y)
{
long n,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;

/* Calculate the number of points */
n = 1;
for (i=0;i<m;i++)
n *= 2;

/* Do the bit reversal */
i2 = n >> 1;
j = 0;
for (i=0;i<n-1;i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}

/* Compute the FFT */
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++) {
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0;j<l1;j++) {
for (i=j;i<n;i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}

/* Scaling for forward transform */
if (dir == 1) {
for (i=0;i<n;i++) {
x[i] /= n;
y[i] /= n;
}
}

return(TRUE);
}

---------------------------------------------------------------------------------
/*****************fft programe*********************/
#include "typedef.h"
#include "math.h"

struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}

void FFT(struct compx *xin,int N)
{
int f,m,nv2,nm1,i,k,j=1,l;
/*int f,m,nv2,nm1,i,k,j=N/2,l;*/
struct compx v,w,t;
nv2=N/2;
f=N;
for(m=1;(f=f/2)!=1;m++){;}
nm1=N-1;

/*變址運算*/
for(i=1;i <=nm1;i++)
{
if(i <j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=nv2;
while(k <j){j=j-k;k=k/2;}
j=j+k;
}

{
int le,lei,ip;
float pi;
for(l=1;l <=m;l++)
{ le=pow(2,l);// 這里用的是L而不是1 !!!!
lei =le/2;
pi=3.14159;
v.real=1.0;
v.imag=0.0;
w.real=cos(pi/lei);
w.imag=-sin(pi/lei);
for(j=1;j <=lei;j++)
{
/*double p=pow(2,m-l)*j;
double ps=2*pi/N*p;
w.real=cos(ps);
w.imag=-sin(ps);*/
for(i=j;i <=N;i=i+le)
{ /* w.real=cos(ps);
w.imag=-sin(ps);*/
ip=i+lei;
t=EE(xin[ip],v);
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
v=EE(v,w);
}
}
}
return;
}

/*****************main programe********************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "typedef.h"

float result[257];
struct compx s[257];
int Num=256;
const float pp=3.14159;

main()
{

int i=1;
for(;i <0x101;i++)
{
s[i].real=sin(pp*i/32);
s[i].imag=0;
}

FFT(s,Num);

for(i=1;i <0x101;i++)
{
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
}

}

-----------------------------------------------------------------------------------
FFT變換 C源代碼

FFT C source code (Simple radix-2)

void fft_float (
unsigned NumSamples,
int InverseTransform,
float *RealIn,
float *ImagIn,
float *RealOut,
float *ImagOut )
{
unsigned NumBits; /* Number of bits needed to store indices */
unsigned i, j, k, n;
unsigned BlockSize, BlockEnd;
double angle_numerator = 2.0 * DDC_PI;
double tr, ti; /* temp real, temp imaginary */
if ( !IsPowerOfTwo(NumSamples) )
{
fprintf (
stderr,
"Error in fft(): NumSamples=%u is not power of two\n",
NumSamples );
exit(1);
}
if ( InverseTransform )
angle_numerator = -angle_numerator;
CHECKPOINTER ( RealIn );
CHECKPOINTER ( RealOut );
CHECKPOINTER ( ImagOut );
NumBits = NumberOfBitsNeeded ( NumSamples );
/*
** Do simultaneous data and bit-reversal ordering into outputs...
*/
for ( i=0; i < NumSamples; i++ )
{
j = ReverseBits ( i, NumBits );
RealOut[j] = RealIn;
ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn;
}
/*
** Do the FFT itself...
*/
BlockEnd = 1;
for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
{
double delta_angle = angle_numerator / (double)BlockSize;
double sm2 = sin ( -2 * delta_angle );
double sm1 = sin ( -delta_angle );
double cm2 = cos ( -2 * delta_angle );
double cm1 = cos ( -delta_angle );
double w = 2 * cm1;
double ar[3], ai[3];
double temp;
for ( i=0; i < NumSamples; i += BlockSize )
{
ar[2] = cm2;
ar[1] = cm1;
ai[2] = sm2;
ai[1] = sm1;
for ( j=i, n=0; n < BlockEnd; j++, n++ )
{
ar[0] = w*ar[1] - ar[2];
ar[2] = ar[1];
ar[1] = ar[0];
ai[0] = w*ai[1] - ai[2];
ai[2] = ai[1];
ai[1] = ai[0];
k = j + BlockEnd;
tr = ar[0]*RealOut[k] - ai[0]*ImagOut[k];
ti = ar[0]*ImagOut[k] + ai[0]*RealOut[k];
RealOut[k] = RealOut[j] - tr;
ImagOut[k] = ImagOut[j] - ti;
RealOut[j] += tr;
ImagOut[j] += ti;
}
}
BlockEnd = BlockSize;
}
/*
** Need to normalize if inverse transform...
*/
if ( InverseTransform )
{
double denom = (double)NumSamples;
for ( i=0; i < NumSamples; i++ )
{
RealOut /= denom;
ImagOut /= denom;
}
}
}

int IsPowerOfTwo ( unsigned x )
{
if ( x < 2 )
return FALSE;
if ( x & (x-1) ) // Thanks to 'byang' for this cute trick!
return FALSE;
return TRUE;
}

unsigned NumberOfBitsNeeded ( unsigned PowerOfTwo )
{
unsigned i;
if ( PowerOfTwo < 2 )
{
fprintf (
stderr,
">>> Error in fftmisc.c: argument %d to NumberOfBitsNeeded is too small.\n",
PowerOfTwo );
exit(1);
}
for ( i=0; ; i++ )
{
if ( PowerOfTwo & (1 << i) )
return i;
}
}

unsigned ReverseBits ( unsigned index, unsigned NumBits )
{
unsigned i, rev;
for ( i=rev=0; i < NumBits; i++ )
{
rev = (rev << 1) | (index & 1);
index >>= 1;
}
return rev;
}

double Index_to_frequency ( unsigned NumSamples, unsigned Index )
{
if ( Index >= NumSamples )
return 0.0;
else if ( Index <= NumSamples/2 )
return (double)Index / (double)NumSamples;
return -(double)(NumSamples-Index) / (double)NumSamples;
}

6. 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變換。

7. 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)來表示