DFT（离散傅里叶变换）

/*******************************************************************************
** 程序名称：离散傅里叶变换(DFT )
** 程序描述：本程序对指定的离散序列进行离散傅里叶变换
** 程序作者：宋元瑞
** 最后修改：2011年4月1日
*******************************************************************************/
#include <stdio.h>
#include <math.h>
#define N	32				//序列长度
#define PI	3.1415926535
typedef double ElementType;	//原始数据序列的数据类型,可以在这里设置
typedef struct				//复数结构体,用于实现傅里叶运算
{
ElementType real,imag;
}complex;
complex dataResult
;		//傅里叶运算的频域结果序列的值(复数)
ElementType dataSource
;	//输入的原始数据序列
ElementType dataFinualResult
; //最终频域序列结果(dataResult
的模,实数)
void FFT_Calculate_OneNode(int k)//计算频域上一个点的DFT值
{
int n = 0;
complex ResultThisNode;
complex part
;
ResultThisNode.real = 0;
ResultThisNode.imag = 0;
for(n=0; n<N; n++)
{
part
.real = cos(2*PI/N*k*n)*dataSource
;//运用欧拉公式把复数拆分成实部和虚部
part
.imag = sin(2*PI/N*k*n)*dataSource
;

ResultThisNode.real += part
.real;
ResultThisNode.imag += part
.imag;
}
dataResult[k].real = ResultThisNode.real;
dataResult[k].imag = ResultThisNode.imag;
}
void FFT_Calculate()//对输入序列全部计算DFT
{
int i = 0;
for(i=0; i<N; i++)
{
FFT_Calculate_OneNode(i);
dataFinualResult[i] = sqrt(dataResult[i].real * dataResult[i].real + dataResult[i].imag * dataResult[i].imag);
}
}
int main(int argc, char *argv[])
{
int i = 0;
for(i=0; i<N; i++)//制造输入序列
{
dataSource[i] = sin(2*PI*i/N);
printf("%lf ",dataSource[i]);
}
FFT_Calculate();//进行DFT计算
printf("/n/n");
for(i=0; i<N; i++)
{printf("%lf ",dataFinualResult[i]);}
return 0;
}