递推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);
}