遞推dft演算法c語言程序

1. 如果這是問題請回答。

如果這是答案請採納。

2. 試根據dft的圖形推導過程，說明dft分析信號的頻譜為什麼只是一種近似的方法

首先，在理解這3個變數之前，你要知道DTFT：

DTFT是離散時間傅里葉變換，用來表達連續的信號的頻譜。

然後理解DFT：

DFT是離散傅里葉變換，針對的是離散的信號和頻譜。DFT是DTFT變化而來，其實就是將連續時間t變成了nT. 為什麼要這樣做呢，因為計算機是在數字環境下工作的，它不可能看見或者處理現實中連續的信號，只能夠進行離散計算，在真實性上盡可能地逼近連續信號。所以DFT是為了我們能夠去用工具分析信號而創造出來的，通常我們直接用DTFT的機會很少。

然後再理解FFT：

記著FFT從本質上來說和DFT沒有任何區別，它只是DFT的一種快速的實現方法而已，比如你要用工具來計算1024個點的DFT來分析一個信號的頻譜，用原來的DFT演算法比起FFT演算法要慢很多，僅此而已。從和硬體的角度看，實現同樣點數的FFT比DFT要快和省程序空間。DSP的書籍都會解釋為什麼FFT實現起來會快一些。

最後理解DCT:

首先，DCT是DFT的一種形式。所謂「餘弦變換」，是在DTFT傅立葉級數展開式中，如果被展開的函數是實偶函數，那麼其傅立葉級數中只包含餘弦項，再將其離散化(DFT)可導出餘弦變換，因此稱之為離散餘弦變換(DCT)。其實DCT屬於DFT的一個子集。DCT用於語音和圖像處理比較多。

希望對你有幫助。

3. 求遞推DFT演算法的C語言程序

風格和地方工會地方化

4. 基於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

5. 用DSP實現FFT演算法變換程序

FFT是DFT的一台計算機與DFT的快速演算法結果
DFT可以說是所有的離散變化的前身是因為類似的變化。
DFT是在時域信號轉換成頻域以一個簡明的物理意義和加工方法的變化。

6. 自己編程實現DFT演算法

(a)
function [Xk]=dft(xn,N)
n=0:1:N-1;
k=0:1:N-1;
WN=exp(-1i*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;

N=256;
n=0:N-1;
xn=cos((5*pi/16)*n);
Xk=dft(xn,N);
subplot(2,1,1)
k=0:1:N-1;
stem(k,abs(Xk));
title('Magnitude of the DFT samples')
xlabel('Frequency index k'); ylabel('Magnitude')
subplot(2,1,2)
stem(k,angle(Xk));
title('Phase of the DFT samples')
xlabel('Frequency index k'); ylabel('Phase')

7. 基數為2的FFT演算法

從上節所述，FFT演算法快速的關鍵在於將原來傅氏矩陣分解為每一行僅含有兩個非零項l與Wⁱ的矩陣的乘積。下面用基數為2，即N=2ⁿ的情形討論矩陣的分解過程.並主要按時間分解的情況討論。

按時間分解的FFT演算法

設N=2ⁿ，n為正整數。考慮輸入序列x₀(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個矩陣，只不過這里的x₁(l)與x₀(l)中的l是用二進制數表示而已。

再寫出(7-2-9)的所有式子，得到

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

將方程組寫成矩陣形式，則有

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

顯然(7-2-15)中的矩陣，正好是(7-1-14)中的方陣，也就是(7-1-12)中被分解出來的第二個矩陣，只不過這里的x₂(l)與x₁(l)是用二進制數表示而已，最後將(7-2-10)中的全部式子寫出，得到

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

將方程組(7-2-16)寫成矩陣形式，則有

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

顯然，(7-2-17)中的方陣正是(7-1-15)中的方陣，也就是(7-1-12)中被分解出來的第1個矩陣，只不過這里的x₃(l)與x₂(l)中的l是用二進制數表示。

由此可見，(7-2-7)中由里往外的三個求和式(7-2-8)、(7-2-9)及(7-2-10)，完全確定了(7-1-12)中三個被分解的矩陣因子。求和得到的最終結果x₃(m₀，m₁，m₂)，與我們所要求的X(m₂，m₁，m₀)正好是逆序的。

到此為止，我們就看到(7-1-11)中的方陣是怎樣被分解成三個方陣因子的。對於N=8，方程(7-2-8)～(7-2-11)就是計算DFT的FFT演算法。為了對FFT演算法有一個直觀的了解並便於編製程序，我們以N=8為例，畫出其流程圖(圖7-2-1)。

根據(7-2-13)，將其中的W⁴用-W⁰代替，畫出從x₀(r)到x₁(r)的流程圖。這一迭代過程用符號#1表示;再根據(7-2-15)，將其中的與W⁴、W⁶分別換成-W⁰與-W²，畫出從x₁(r)到x₂(r)的流程圖，這一迭代過程記為#2;最後，根據(7-2-17)，將其中的W⁴、W⁶、W⁵、W⁷分別換成-W⁰、-W²、-W¹、-W³，畫出流程圖7-2-3合並圖7-2-1～7-2-3，就得到從x₀(r)到x₃(r)的完整流程圖7-2-4。

圖7-2-1 第一次遞推

圖7-2-2 第二次遞推

在圖7-2-5中，畫出N=16=2⁴的FFT演算法流程圖:

根據從x₀到譜X的FFT演算法流程圖7-2-4與圖7-2-5，我們總結出如下幾點:

(1)從x₀到終值x_r的最大迭代次數r，由r=log₂N確定。

例如，N=8，最大迭代次數r=3;N=16，最大迭代次數r=4。

(2)在第r次迭代中，要乘的W因子為

圖7-2-3 第三次遞推

例如，N=8，在第一次迭代中，要乘因子W⁰;在第二次迭代中要乘因子W⁰，W¹，W²，W³。

(3)在第r次迭代中，包含2^r-1個組，每個組

包含 。例如N=8，第一次迭代r=1，有

一個組，每組包含8個x(s);在第二次迭代中包含2個組，每組包含4個x(s);第三次迭代中包含4個組，每組2個x(s)。

圖7-2-4 x₀(r)到x₃(r)的完整流程圖

(4)在第r次迭代中，各組包含的W因子各不相同，且每一組僅包含一種類型的因子 ，此因子在組中一半數為正，另一半數為負。例如N=8，第二次迭代中，第一組包含因子W⁰，且在該組中一半數為正，另一半數為負;第二組包含因子W²，在該組中也是一半數為正，另一半數為負。

(5)在第r次迭代中，各組包含的W因子除正負號外類型均相同。所以只須確定每組中第一個因子，之後按半數反號，即得到所求W因子。具體做法是，在每組第一個因子W^SN2r對應的x_r(k)中，將k表示成n位的二進制數，n=log₂N，並把這個二進制數右移n-r位，左邊空出的地方添零補足n位，之後再將此n位二進制數逆位，即得到所求W因子的指數。例如，N=8，迭代#2包含兩組，每組包含4個x₂(k)，第二組第一個因子W對應於x₂(4)。將4表示成3位的二進制數為100，把它右移1位成10，右邊添零補成3位為010，逆位仍為010，故所求因子為W²，第2組全部W因子為W²，W²，-W²，-W²。又如，N=16，迭代#3中包含4個組，每組包含4個x₃(k)，第4組第1個因子W對應於x₃(12)。將12表示成4位的二進制數為1100。右移1位變成110，將左邊空處添零補成4位為0110，逆位仍為0110，故所求因子為W⁶，從而第4組全部W因子為W⁶，W⁶，-W⁶與-W⁶。

圖7-2-5 N=16=2⁴的FFT演算法流程圖

(6)如果已知N=2的FFT演算法，容易從它求得n=2的FFT演算法。具體作法是，在n=2^n-1的FFT演算法中，將所有x_r(l)的個數加倍，所有W的個數及其乘冪加倍，就得N=2ⁿ中前n-1次迭代的全部結果。之後，將新得到的第n-1次迭代中乘冪相同的W個數減半，就是第n次迭代中前2ⁿ/2個W，將這些W的乘冪依次加1，就得到後2ⁿ/2個W。例如N=16的前3次迭代，都是N=8的三次迭代中所有x_r(l)的個數加倍，W的個數及其乘冪加倍的結果。再將N=16的第三次迭代中乘冪相同的W個數減半，就是第4次迭代中前8個W。

(7)在第r-1次迭代中的x_r-1(i)與x_r-1(j)僅用於計算r次迭代中的x_r(i)與x_r(j)，而不會用於計算任何其它的x_r(k)與x_r(l)。例如N=16的第二次迭代中的x₂(0)與x₂(2)，只用於計算第三次迭代中的x₃(0)與x₃(2);第三次迭代中的x₃(8)與x₃(9)也只用於計算第四次迭代中的x₄(8)與x₄(9)。因此，我們可以把第r次迭代中的x_r(i)與x_r(j)，存放到r-1次迭代x_r-1(i)與x_r-1(j)所佔用的原存儲單元中去，這樣，所需要的計算機存儲容量就只限於原數據序列占據的單元數。如果是復序列的話，存儲單元要加倍。

(8)上述FFT演算法也可用於計算逆離散傅氏變換(IDFT)(圖7-2-6)，只不過在計算時要將上述FFT演算法中的W因子用其共軛W^*代替，並將最後迭代計算的結果統統乘以1/N.

圖7-2-6 N=8的逆離散富氏變換流程圖

8. 求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 中

//引用任何所需的附加頭文件，而不是在此文件中引用

9. 設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;

(9)遞推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為基底的混合基演算法。

10. 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);
}