⑴ 基數為2的FFT演算法
從上節所述，FFT演算法快速的關鍵在於將原來傅氏矩陣分解為每一行僅含有兩個非零項l與Wi的矩陣的乘積。下面用基數為2，即N=2n的情形討論矩陣的分解過程.並主要按時間分解的情況討論。
按時間分解的FFT演算法
設N=2n，n為正整數。考慮輸入序列x0(l)的DFT
物探數字信號分析與處理技術
將l與m用二進製表示
物探數字信號分析與處理技術
將(7-2-2)代入(7-2-1)中，得到
物探數字信號分析與處理技術
為了說明問題，我們取N=8，於是從(7-2-2)得到
物探數字信號分析與處理技術
從(7-2-4)和(7-2-3)得到
物探數字信號分析與處理技術
將(7-2-5)中的W的指數按時間l進行分解，有
物探數字信號分析與處理技術
因為
物探數字信號分析與處理技術
故從(7-2-6)得到
物探數字信號分析與處理技術
將上式代入(7-2-5)中得到
物探數字信號分析與處理技術
物探數字信號分析與處理技術
我們在公式(7-2-7)中由里往外求和，並置
物探數字信號分析與處理技術
於是得到
物探數字信號分析與處理技術
首先寫出(7-2-8)的所有式子
物探數字信號分析與處理技術
將方程組(7-2-12)寫成矩陣形式，得到
物探數字信號分析與處理技術
我們看到(7-2-13)中的方陣，正好是(7-1-13)中的方陣，也就是(7-1-12)中被分解出來的第3個矩陣，只不過這里的x1(l)與x0(l)中的l是用二進制數表示而已。
再寫出(7-2-9)的所有式子，得到
物探數字信號分析與處理技術
將方程組寫成矩陣形式，則有
物探數字信號分析與處理技術
顯然(7-2-15)中的矩陣，正好是(7-1-14)中的方陣，也就是(7-1-12)中被分解出來的第二個矩陣，只不過這里的x2(l)與x1(l)是用二進制數表示而已，最後將(7-2-10)中的全部式子寫出，得到
物探數字信號分析與處理技術
將方程組(7-2-16)寫成矩陣形式，則有
物探數字信號分析與處理技術
顯然，(7-2-17)中的方陣正是(7-1-15)中的方陣，也就是(7-1-12)中被分解出來的第1個矩陣，只不過這里的x3(l)與x2(l)中的l是用二進制數表示。
由此可見，(7-2-7)中由里往外的三個求和式(7-2-8)、(7-2-9)及(7-2-10)，完全確定了(7-1-12)中三個被分解的矩陣因子。求和得到的最終結果x3(m0，m1，m2)，與我們所要求的X(m2，m1，m0)正好是逆序的。
到此為止，我們就看到(7-1-11)中的方陣是怎樣被分解成三個方陣因子的。對於N=8，方程(7-2-8)～(7-2-11)就是計算DFT的FFT演算法。為了對FFT演算法有一個直觀的了解並便於編製程序，我們以N=8為例，畫出其流程圖(圖7-2-1)。
根據(7-2-13)，將其中的W4用-W0代替，畫出從x0(r)到x1(r)的流程圖。這一迭代過程用符號#1表示;再根據(7-2-15)，將其中的與W4、W6分別換成-W0與-W2，畫出從x1(r)到x2(r)的流程圖，這一迭代過程記為#2;最後，根據(7-2-17)，將其中的W4、W6、W5、W7分別換成-W0、-W2、-W1、-W3，畫出流程圖7-2-3合並圖7-2-1～7-2-3，就得到從x0(r)到x3(r)的完整流程圖7-2-4。
圖7-2-1 第一次遞推
圖7-2-2 第二次遞推
在圖7-2-5中，畫出N=16=24的FFT演算法流程圖:
根據從x0到譜X的FFT演算法流程圖7-2-4與圖7-2-5，我們總結出如下幾點:
(1)從x0到終值xr的最大迭代次數r，由r=log2N確定。
例如，N=8，最大迭代次數r=3;N=16，最大迭代次數r=4。
(2)在第r次迭代中，要乘的W因子為
圖7-2-3 第三次遞推
例如，N=8，在第一次迭代中，要乘因子W0;在第二次迭代中要乘因子W0，W1，W2，W3。
(3)在第r次迭代中，包含2r-1個組，每個組
包含 。例如N=8，第一次迭代r=1，有
一個組，每組包含8個x(s);在第二次迭代中包含2個組，每組包含4個x(s);第三次迭代中包含4個組，每組2個x(s)。
圖7-2-4 x0(r)到x3(r)的完整流程圖
(4)在第r次迭代中，各組包含的W因子各不相同，且每一組僅包含一種類型的因子 ，此因子在組中一半數為正，另一半數為負。例如N=8，第二次迭代中，第一組包含因子W0，且在該組中一半數為正，另一半數為負;第二組包含因子W2，在該組中也是一半數為正，另一半數為負。
(5)在第r次迭代中，各組包含的W因子除正負號外類型均相同。所以只須確定每組中第一個因子，之後按半數反號，即得到所求W因子。具體做法是，在每組第一個因子WSN2r對應的xr(k)中，將k表示成n位的二進制數，n=log2N，並把這個二進制數右移n-r位，左邊空出的地方添零補足n位，之後再將此n位二進制數逆位，即得到所求W因子的指數。例如，N=8，迭代#2包含兩組，每組包含4個x2(k)，第二組第一個因子W對應於x2(4)。將4表示成3位的二進制數為100，把它右移1位成10，右邊添零補成3位為010，逆位仍為010，故所求因子為W2，第2組全部W因子為W2，W2，-W2，-W2。又如，N=16，迭代#3中包含4個組，每組包含4個x3(k)，第4組第1個因子W對應於x3(12)。將12表示成4位的二進制數為1100。右移1位變成110，將左邊空處添零補成4位為0110，逆位仍為0110，故所求因子為W6，從而第4組全部W因子為W6，W6，-W6與-W6。
圖7-2-5 N=16=24的FFT演算法流程圖
(6)如果已知N=2的FFT演算法，容易從它求得n=2的FFT演算法。具體作法是，在n=2n-1的FFT演算法中，將所有xr(l)的個數加倍，所有W的個數及其乘冪加倍，就得N=2n中前n-1次迭代的全部結果。之後，將新得到的第n-1次迭代中乘冪相同的W個數減半，就是第n次迭代中前2n/2個W，將這些W的乘冪依次加1，就得到後2n/2個W。例如N=16的前3次迭代，都是N=8的三次迭代中所有xr(l)的個數加倍，W的個數及其乘冪加倍的結果。再將N=16的第三次迭代中乘冪相同的W個數減半，就是第4次迭代中前8個W。
(7)在第r-1次迭代中的xr-1(i)與xr-1(j)僅用於計算r次迭代中的xr(i)與xr(j)，而不會用於計算任何其它的xr(k)與xr(l)。例如N=16的第二次迭代中的x2(0)與x2(2)，只用於計算第三次迭代中的x3(0)與x3(2);第三次迭代中的x3(8)與x3(9)也只用於計算第四次迭代中的x4(8)與x4(9)。因此，我們可以把第r次迭代中的xr(i)與xr(j)，存放到r-1次迭代xr-1(i)與xr-1(j)所佔用的原存儲單元中去，這樣，所需要的計算機存儲容量就只限於原數據序列占據的單元數。如果是復序列的話，存儲單元要加倍。
(8)上述FFT演算法也可用於計算逆離散傅氏變換(IDFT)(圖7-2-6)，只不過在計算時要將上述FFT演算法中的W因子用其共軛W*代替，並將最後迭代計算的結果統統乘以1/N.
圖7-2-6 N=8的逆離散富氏變換流程圖
⑵ 基於FFT的演算法優化 要c語言完整程序（利用旋轉因子的性質），有的請留言，答謝！！！(有核心代碼，望指教
實現(C描述)
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
//#include "complex.h"
// --------------------------------------------------------------------------
#define N 8 //64
#define M 3 //6 //2^m=N
#define PI 3.1415926
// --------------------------------------------------------------------------
float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};
float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};
float x_i[N]; //N=8
/*
float twiddle[N/2] = {1, 0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,
0.7075, 0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991,
0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,
-0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951}; //N=64
float x_r[N]={1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,};
float x_i[N];
*/
FILE *fp;
// ----------------------------------- func -----------------------------------
/**
* 初始化輸出虛部
*/
static void fft_init( void )
{
int i;
for(i=0; i<N; i++) x_i[i] = 0.0;
}
/**
* 反轉演算法.將時域信號重新排序.
* 這個演算法有改進的空間
*/
static void bitrev( void )
{
int p=1, q, i;
int bit_rev[ N ]; //
float xx_r[ N ]; //
bit_rev[ 0 ] = 0;
while( p < N )
{
for(q=0; q<p; q++)
{
bit_rev[ q ] = bit_rev[ q ] * 2;
bit_rev[ q + p ] = bit_rev[ q ] + 1;
}
p *= 2;
}
for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];
for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ];
}
/* ------------ add by sshc625 ------------ */
static void bitrev2( void )
{
return ;
}
/* */
void display( void )
{
printf("\n\n");
int i;
for(i=0; i<N; i++)
printf("%f\t%f\n", x_r[i], x_i[i]);
}
/**
*
*/
void fft1( void )
{ fp = fopen("log1.txt", "a+");
int L, i, b, j, p, k, tx1, tx2;
float TR, TI, temp; // 臨時變數
float tw1, tw2;
/* 深M. 對層進行循環. L為當前層, 總層數為M. */
for(L=1; L<=M; L++)
{
fprintf(fp,"----------Layer=%d----------\n", L);
/* b的意義非常重大,b表示當前層的顆粒具有的輸入樣本點數 */
b = 1;
i = L - 1;
while(i > 0)
{
b *= 2;
i--;
}
// -------------- 是否外層對顆粒循環, 內層對樣本點循環邏輯性更強一些呢! --------------
/*
* outter對參與DFT的樣本點進行循環
* L=1, 循環了1次(4個顆粒, 每個顆粒2個樣本點)
* L=2, 循環了2次(2個顆粒, 每個顆粒4個樣本點)
* L=3, 循環了4次(1個顆粒, 每個顆粒8個樣本點)
*/
for(j=0; j<b; j++)
{
/* 求旋轉因子tw1 */
p = 1;
i = M - L; // M是為總層數, L為當前層.
while(i > 0)
{
p = p*2;
i--;
}
p = p * j;
tx1 = p % N;
tx2 = tx1 + 3*N/4;
tx2 = tx2 % N;
// tw1是cos部分, 實部; tw2是sin部分, 虛數部分.
tw1 = ( tx1>=N/2)? -twiddle[tx1-N/2] : twiddle[ tx1 ];
tw2 = ( tx2>=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];
/*
* inner對顆粒進行循環
* L=1, 循環了4次(4個顆粒, 每個顆粒2個輸入)
* L=2, 循環了2次(2個顆粒, 每個顆粒4個輸入)
* L=3, 循環了1次(1個顆粒, 每個顆粒8個輸入)
*/
for(k=j; k<N; k=k+2*b)
{
TR = x_r[k]; // TR就是A, x_r[k+b]就是B.
TI = x_i[k];
temp = x_r[k+b];
/*
* 如果復習一下 (a+j*b)(c+j*d)兩個復數相乘後的實部虛部分別是什麼
* 就能理解為什麼會如下運算了, 只有在L=1時候輸入才是實數, 之後層的
* 輸入都是復數, 為了讓所有的層的輸入都是復數, 我們只好讓L=1時候的
* 輸入虛部為0
* x_i[k+b]*tw2是兩個虛數相乘
*/
fprintf(fp, "tw1=%f, tw2=%f\n", tw1, tw2);
x_r[k] = TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;
x_i[k] = TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;
x_r[k+b] = TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;
x_i[k+b] = TI + temp*tw2 - x_i[k+b]*tw1;
fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f\n", k, x_r[k], x_i[k]);
fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f\n", k+b, x_r[k+b], x_i[k+b]);
} //
} //
} //
}
/**
* ------------ add by sshc625 ------------
* 該實現的流程為
* for( Layer )
* for( Granule )
* for( Sample )
*
*
*
*
*/
void fft2( void )
{ fp = fopen("log2.txt", "a+");
int cur_layer, gr_num, i, k, p;
float tmp_real, tmp_imag, temp; // 臨時變數, 記錄實部
float tw1, tw2;// 旋轉因子,tw1為旋轉因子的實部cos部分, tw2為旋轉因子的虛部sin部分.
int step; // 步進
int sample_num; // 顆粒的樣本總數(各層不同, 因為各層顆粒的輸入不同)
/* 對層循環 */
for(cur_layer=1; cur_layer<=M; cur_layer++)
{
/* 求當前層擁有多少個顆粒(gr_num) */
gr_num = 1;
i = M - cur_layer;
while(i > 0)
{
i--;
gr_num *= 2;
}
/* 每個顆粒的輸入樣本數N' */
sample_num = (int)pow(2, cur_layer);
/* 步進. 步進是N'/2 */
step = sample_num/2;
/* */
k = 0;
/* 對顆粒進行循環 */
for(i=0; i<gr_num; i++)
{
/*
* 對樣本點進行循環, 注意上限和步進
*/
for(p=0; p<sample_num/2; p++)
{
// 旋轉因子, 需要優化...
tw1 = cos(2*PI*p/pow(2, cur_layer));
tw2 = -sin(2*PI*p/pow(2, cur_layer));
tmp_real = x_r[k+p];
tmp_imag = x_i[k+p];
temp = x_r[k+p+step];
/*(tw1+jtw2)(x_r[k]+jx_i[k])
*
* real : tw1*x_r[k] - tw2*x_i[k]
* imag : tw1*x_i[k] + tw2*x_r[k]
* 我想不抽象出一個
* typedef struct {
* double real; // 實部
* double imag; // 虛部
* } complex; 以及針對complex的操作
* 來簡化復數運算是否是因為效率上的考慮!
*/
/* 蝶形演算法 */
x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );
x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );
/* X[k] = A(k)+WB(k)
* X[k+N/2] = A(k)-WB(k) 的性質可以優化這里*/
// 旋轉因子, 需要優化...
tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));
tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));
x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );
x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );
printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);
printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);
}
/* 開跳!:) */
k += 2*step;
}
}
}
/*
* 後記:
* 究竟是顆粒在外層循環還是樣本輸入在外層, 好象也差不多, 復雜度完全一樣.
* 但以我資質愚鈍花費了不少時間才弄明白這數十行代碼.
* 從中我發現一個於我非常有幫助的教訓, 很久以前我寫過一部分演算法, 其中絕大多數都是遞歸.
* 將數據量減少, 減少再減少, 用歸納的方式來找出數據量加大代碼的規律
* 比如FFT
* 1. 先寫死LayerI的代碼; 然後再把LayerI的輸出作為LayerII的輸入, 又寫死代碼; ......
* 大約3層就可以統計出規律來. 這和遞歸也是一樣, 先寫死一兩層, 自然就出來了!
* 2. 有的功能可以寫偽代碼, 不急於求出結果, 降低復雜性, 把邏輯結果定出來後再添加.
* 比如旋轉因子就可以寫死, 就寫1.0. 流程出來後再寫旋轉因子.
* 寥寥數語, 我可真是流了不少汗! Happy!
*/
void dft( void )
{
int i, n, k, tx1, tx2;
float tw1,tw2;
float xx_r[N],xx_i[N];
/*
* clear any data in Real and Imaginary result arrays prior to DFT
*/
for(k=0; k<=N-1; k++)
xx_r[k] = xx_i[k] = x_i[k] = 0.0;
// caculate the DFT
for(k=0; k<=(N-1); k++)
{
for(n=0; n<=(N-1); n++)
{
tx1 = (n*k);
tx2 = tx1+(3*N)/4;
tx1 = tx1%(N);
tx2 = tx2%(N);
if(tx1 >= (N/2))
tw1 = -twiddle[tx1-(N/2)];
else
tw1 = twiddle[tx1];
if(tx2 >= (N/2))
tw2 = -twiddle[tx2-(N/2)];
else
tw2 = twiddle[tx2];
xx_r[k] = xx_r[k]+x_r[n]*tw1;
xx_i[k] = xx_i[k]+x_r[n]*tw2;
}
xx_i[k] = -xx_i[k];
}
// display
for(i=0; i<N; i++)
printf("%f\t%f\n", xx_r[i], xx_i[i]);
}
// ---------------------------------------------------------------------------
int main( void )
{
fft_init( );
bitrev( );
// bitrev2( );
//fft1( );
fft2( );
display( );
system( "pause" );
// dft();
return 1;
}
本文來自CSDN博客，轉載請標明出處：http://blog.csdn.net/sshcx/archive/2007/06/14/1651616.aspx
⑶ 設x(n)={1，0.5，0，0.5，1，1，0.5，0)，用FFT演算法求x(n)的DFT。FFT演算法任選，畫出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行。
例：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 1000
/*定義復數類型*/
typedef struct{
double real;
double img;
}complex;
complex x[N], *W; /*輸入序列,變換核*/
int size_x=0;/*輸入序列的大小，在本程序中僅限2的次冪*/
double PI;/*圓周率*/
void fft();/*快速傅里葉變換*/
void initW(); /*初始化變換核*/
void change(); /*變址*/
void add(complex ,complex ,complex *); /*復數加法*/
void mul(complex ,complex ,complex *); /*復數乘法*/
void sub(complex ,complex ,complex *); /*復數減法*/
void output();
int main(){
int i;/*輸出結果*/
system("cls");
PI=atan(1)*4;
printf("Please input the size of x: ");
scanf("%d",&size_x);
printf("Please input the data in x[N]: ");
for(i=0;i<size_x;i++)
scanf("%lf%lf",&x[i].real,&x[i].img);
initW();
fft();
output();
return 0;
}
/*快速傅里葉變換*/
void fft(){
int i=0,j=0,k=0,l=0;
complex up,down,proct;
change();
for(i=0;i< log(size_x)/log(2) ;i++){ /*一級蝶形運算*/
l=1<<i;
(3)遞推dft演算法c語言程序擴展閱讀：
FFT演算法很多，根據實現運算過程是否有指數因子WN可分為有、無指數因子的兩類演算法。
經典庫利-圖基演算法 當輸入序列的長度N不是素數(素數只能被1而它本身整除）而是可以高度分解的復合數，即N=N1N2N3…Nr時，若N1=N2=…=Nr=2，N=2則N點DFT的計算可分解為N=2×N/2，即兩個N/2點DFT計算的組合，而N/2點DFT的計算又可分解為N/2=2×N/4，即兩個N/4點DFT計算的組合。
依此類推，使DFT的計算形成有規則的模式，故稱之為以2為基底的FFT演算法。同理，當N=4時，則稱之為以4為基底的FFT演算法。當N=N1·N2時，稱為以N1和N2為基底的混合基演算法。
⑷ 求 原始DFT 演算法 代碼。有的麻煩發到我郵箱[email protected]
加入多個PushButton，如何設置始終有個獲得焦點（這個比較好實現）
但是，如下情況的時候要怎麼處理：一個界面中左邊是多個PushButton，右邊也是可獲得焦點的組
⑸ C語言程序，離散傅里葉變換，調用函數聲明處總有錯誤
int DFT(int dir,int m,double *x1,double *y1)
{
long i,k;
double arg;
double cosarg,sinarg;
double *x2=NULL,*y2=NULL;
x2=malloc(m*sizeof(double));
y2=malloc(m*sizeof(double));
if(x2==NULL||y2==NULL)return(FALSE);
for(i=0;i<m;i++)
{
x2[i]=0;
y2[i]=0;
arg=-dir*2.0*3.141592654*(double)i/(double)m;
for(k=0;k<m;k++)
{
cosarg=cos(k*arg);
sinarg=sin(k*arg);
x2[i]+=(x1[k]*cosarg-y1[k]*sinarg);
y2[i]+=(x1[k]*sinarg+y1[k]*cosarg);
}
}
/*Copythedataback*/
if(dir==1)
{
for(i=0;i<m;i++)
{
x1[i]=x2[i]/(double)m;
y1[i]=y2[i]/(double)m;
}
}
else
{
for(i=0;i<m;i++)
{
x1[i]=x2[i];
y1[i]=y2[i];
}
}
free(x2);
free(y2);
return(TRUE);
}
⑹ 求遞推DFT演算法的C語言程序
風格和地方工會地方化
⑺ 我做「基於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:
⑻ DFT變換演算法的全稱
離散傅里葉變換(Discrete Fourier Transform,縮寫為DFT),是傅里葉變換在時域和頻域上都呈離散的形式,將信號的時域采樣變換為其DTFT的頻域描述了離散傅里葉變換的推導過程,是計算機相關演算法編寫的依據,;並且可因此推出快速傅里葉變換演算法
⑼ 求FFT的c語言程序
快速傅里葉變換 要用C++ 才行吧 你可以用MATLAB來實現更方便點啊
此FFT 是用VC6.0編寫，由FFT.CPP；STDAFX.H和STDAFX.CPP三個文件組成，編譯成功。程序可以用文件輸入和輸出為文件。文件格式為TXT文件。測試結果如下：
輸入文件：8.TXT 或手動輸入
8 //N
1
2
3
4
5
6
7
8
輸出結果為：或保存為TXT文件。（8OUT.TXT）
8
(36,0)
(-4,9.65685)
(-4,4)
(-4,1.65685)
(-4,0)
(-4,-1.65685)
(-4,-4)
(-4,-9.65685)
下面為FFT.CPP文件：
// FFT.cpp : 定義控制台應用程序的入口點。
#include "stdafx.h"
#include <iostream>
#include <complex>
#include <bitset>
#include <vector>
#include <conio.h>
#include <string>
#include <fstream>
using namespace std;
bool inputData(unsigned long &, vector<complex<double> >&); //手工輸入數據
void FFT(unsigned long &, vector<complex<double> >&); //FFT變換
void display(unsigned long &, vector<complex<double> >&); //顯示結果
bool readDataFromFile(unsigned long &, vector<complex<double> >&); //從文件中讀取數據
bool saveResultToFile(unsigned long &, vector<complex<double> >&); //保存結果至文件中
const double PI = 3.1415926;
int _tmain(int argc, _TCHAR* argv[])
{
vector<complex<double> > vecList; //有限長序列
unsigned long ulN = 0; //N
char chChoose = ' '; //功能選擇
//功能循環
while(chChoose != 'Q' && chChoose != 'q')
{
//顯示選擇項
cout << "\nPlease chose a function" << endl;
cout << "\t1.Input data manually, press 'M':" << endl;
cout << "\t2.Read data from file, press 'F':" << endl;
cout << "\t3.Quit, press 'Q'" << endl;
cout << "Please chose:";
//輸入選擇
chChoose = getch();
//判斷
switch(chChoose)
{
case 'm': //手工輸入數據
case 'M':
if(inputData(ulN, vecList))
{
FFT(ulN, vecList);
display(ulN, vecList);
saveResultToFile(ulN, vecList);
}
break;
case 'f': //從文檔讀取數據
case 'F':
if(readDataFromFile(ulN, vecList))
{
FFT(ulN, vecList);
display(ulN, vecList);
saveResultToFile(ulN, vecList);
}
break;
}
}
return 0;
}
bool Is2Power(unsigned long ul) //判斷是否是2的整數次冪
{
if(ul < 2)
return false;
while( ul > 1 )
{
if( ul % 2 )
return false;
ul /= 2;
}
return true;
}
bool inputData(unsigned long & ulN, vector<complex<double> >& vecList)
{
//題目
cout<< "\n\n\n==============================Input Data===============================" << endl;
//輸入N
cout<< "\nInput N:";
cin>>ulN;
if(!Is2Power(ulN)) //驗證N的有效性
{
cout<< "N is invalid (N must like 2, 4, 8, .....), please retry." << endl;
return false;
}
//輸入各元素
vecList.clear(); //清空原有序列
complex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
cout << "Input x(" << i << "):";
cin >> c;
vecList.push_back(c);
}
return true;
}
bool readDataFromFile(unsigned long & ulN, vector<complex<double> >& vecList) //從文件中讀取數據
{
//題目
cout<< "\n\n\n===============Read Data From File==============" << endl;
//輸入文件名
string strfilename;
cout << "Input filename:" ;
cin >> strfilename;
//打開文件
cout << "open file " << strfilename << "......." <<endl;
ifstream loadfile;
loadfile.open(strfilename.c_str());
if(!loadfile)
{
cout << "\tfailed" << endl;
return false;
}
else
{
cout << "\tsucceed" << endl;
}
vecList.clear();
//讀取N
loadfile >> ulN;
if(!loadfile)
{
cout << "can't get N" << endl;
return false;
}
else
{
cout << "N = " << ulN << endl;
}
//讀取元素
complex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
loadfile >> c;
if(!loadfile)
{
cout << "can't get enough infomation" << endl;
return false;
}
else
cout << "x(" << i << ") = " << c << endl;
vecList.push_back(c);
}
//關閉文件
loadfile.close();
return true;
}
bool saveResultToFile(unsigned long & ulN, vector<complex<double> >& vecList) //保存結果至文件中
{
//詢問是否需要將結果保存至文件
char chChoose = ' ';
cout << "Do you want to save the result to file? (y/n):";
chChoose = _getch();
if(chChoose != 'y' && chChoose != 'Y')
{
return true;
}
//輸入文件名
string strfilename;
cout << "\nInput file name:" ;
cin >> strfilename;
cout << "Save result to file " << strfilename << "......" << endl;
//打開文件
ofstream savefile(strfilename.c_str());
if(!savefile)
{
cout << "can't open file" << endl;
return false;
}
//寫入N
savefile << ulN << endl;
//寫入元素
for(vector<complex<double> >::iterator i = vecList.begin(); i < vecList.end(); i++)
{
savefile << *i << endl;
}
//寫入完畢
cout << "save succeed." << endl;
//關閉文件
savefile.close();
return true;
}
void FFT(unsigned long & ulN, vector<complex<double> >& vecList)
{
//得到冪數
unsigned long ulPower = 0; //冪數
unsigned long ulN1 = ulN - 1;
while(ulN1 > 0)
{
ulPower++;
ulN1 /= 2;
}
//反序
bitset<sizeof(unsigned long) * 8> bsIndex; //二進制容器
unsigned long ulIndex; //反轉後的序號
unsigned long ulK;
for(unsigned long p = 0; p < ulN; p++)
{
ulIndex = 0;
ulK = 1;
bsIndex = bitset<sizeof(unsigned long) * 8>(p);
for(unsigned long j = 0; j < ulPower; j++)
{
ulIndex += bsIndex.test(ulPower - j - 1) ? ulK : 0;
ulK *= 2;
}
if(ulIndex > p)
{
complex<double> c = vecList[p];
vecList[p] = vecList[ulIndex];
vecList[ulIndex] = c;
}
}
//計算旋轉因子
vector<complex<double> > vecW;
for(unsigned long i = 0; i < ulN / 2; i++)
{
vecW.push_back(complex<double>(cos(2 * i * PI / ulN) , -1 * sin(2 * i * PI / ulN)));
}
for(unsigned long m = 0; m < ulN / 2; m++)
{
cout<< "\nvW[" << m << "]=" << vecW[m];
}
//計算FFT
unsigned long ulGroupLength = 1; //段的長度
unsigned long ulHalfLength = 0; //段長度的一半
unsigned long ulGroupCount = 0; //段的數量
complex<double> cw; //WH(x)
complex<double> c1; //G(x) + WH(x)
complex<double> c2; //G(x) - WH(x)
for(unsigned long b = 0; b < ulPower; b++)
{
ulHalfLength = ulGroupLength;
ulGroupLength *= 2;
for(unsigned long j = 0; j < ulN; j += ulGroupLength)
{
for(unsigned long k = 0; k < ulHalfLength; k++)
{
cw = vecW[k * ulN / ulGroupLength] * vecList[j + k + ulHalfLength];
c1 = vecList[j + k] + cw;
c2 = vecList[j + k] - cw;
vecList[j + k] = c1;
vecList[j + k + ulHalfLength] = c2;
}
}
}
}
void display(unsigned long & ulN, vector<complex<double> >& vecList)
{
cout << "\n\n===========================Display The Result=========================" << endl;
for(unsigned long d = 0; d < ulN;d++)
{
cout << "X(" << d << ")\t\t\t = " << vecList[d] << endl;
}
}
下面為STDAFX.H文件：
// stdafx.h : 標准系統包含文件的包含文件，
// 或是常用但不常更改的項目特定的包含文件
#pragma once
#include <iostream>
#include <tchar.h>
// TODO: 在此處引用程序要求的附加頭文件
下面為STDAFX.CPP文件：
// stdafx.cpp : 只包括標准包含文件的源文件
// FFT.pch 將成為預編譯頭
// stdafx.obj 將包含預編譯類型信息
#include "stdafx.h"
// TODO: 在 STDAFX.H 中
//引用任何所需的附加頭文件，而不是在此文件中引用
⑽ 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;
}