C语言编写DFT计算程序， 并绘制幅度谱

一、什么是DFT？
离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作其周期延拓的变换。在实际应用中通常采用快速傅里叶变换计算DFT

二、DFT的公式
离散傅里叶变化定义为
$$X(k)=\sum _{n=0}^{N-1}x(n)W_N^{nk},k=0,1,\dots ,N-1$$

其中，$$W_N^{nk}$$为旋转因子，特此说明，其计算公式为$$W_N^{nk}=e^{-j\frac{2\pi }{N}nk},\mathrm{N}=0,1_{\mu },\dots ,\mathrm{N}-1$$
1.利用欧拉公式$$e^{j\theta }=\cos \theta +j\sin \theta$$得：$$X(k)=\sum _{n=0}^{N-1}x(n)\cos \left(\frac{2\pi }{N}kn\right)-jx(n)\sin \left(\frac{2\pi }{N}kn\right)$$
2.幅值得计算
$$\mathrm{asin}(x)+b\cos (x)=\sqrt{a^2+b^2}\sin (x+\mathrm{arctanb}/a)$$

三，DFT源程序
抽样点为：64
原始序列为：
0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)

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

#define pi 3.14159
#define N 64

typedef struct
{
	double real;
	double imag;
}complex;

void Wn(int ,int ,int ,complex *);//定义旋转因子
void c_jiafa(complex ,complex ,complex *); //复数加法运算 
void c_chengfa(complex ,complex ,complex *); //复数乘法运算

int main()
{
	double y[N];
	complex A[N],a[N],f[N]; //a[N]为输入的信号，A[N]是DFT后的序列 
	for(int i=0;i<N;i++) //提供初始信号 
	{
		a[i].real=0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i);//0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)
		a[i].imag=0;
	}
	
	
	for(int k=0;k<N;k++)
	{
		//int n;
		complex sum={0,0};
		for(int n=0;n<N;n++)
		{
			complex wn,t;
			Wn(N,n,k,&wn);
			c_chengfa(a[n],wn,&t);
			c_jiafa(sum,t,&sum);
		}
		A[k].real=sum.real;
		A[k].imag=sum.imag;
	} 
	
	for(int i=0;i<N;i++)//幅值计算
	{
		y[i]=sqrt(A[i].real*A[i].real+A[i].imag*A[i].imag);
		printf("%d\t%lf\n",i,y[i]);
	}
	system("pause");
	return 0;
}

void Wn(int n1,int i,int k,complex *Wn) //定义旋转因子 
{
	//用欧拉公式分成实部和虚部 
	Wn->real=cos(2*pi*i*k/n1);
	Wn->imag=-sin(2*pi*i*k/n1);
}

void c_jiafa(complex a,complex b,complex *c) //复数加法运算 
{
	c->real=a.real+b.real; 
	c->imag=a.imag+b.imag;
}

void c_chengfa(complex a,complex b,complex *c) //复数乘法运算
{
	c->real=a.real*b.real-a.imag*b.imag;
	c->imag=a.real*b.imag+a.imag*b.real;
}

四，程序编译，产生dat文件

五，gnuplot绘图

结果为：