C语言有限长序列的DFT实现及gnuplot作图

离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是傅里叶变换在时域和频域上都呈离散的形式。在形式上，时域和频域上的序列都是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。对于有限长的离散信号作DFT，可以把它看作是离散周期信号的一个周期，应进行周期延拓之后再进行计算。

设有限长序列的长度为N，那么其离散傅里叶变换定义为：

我们用C语言实现DFT计算，实际就是求f(n)的值，也就是求信号的幅值。首先我们使用到指数的欧拉公式：

代入DFT的计算式可以得到：

要计算最后的幅值，需要将信号分为实部和虚部，再由正余弦的合并公式：

我们可以知道，f(n)的实部和虚部分别表示为：

那么幅值的计算式为：

程序代码：

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

#define PI 3.14159

struct data
{
	double real,imag;
};

int main()
{
	int N=16;                  //序列长度

	double Squence[100];       //序列
	double Amplitude[100];     //幅值

	struct data one[16];       //单个数据
	struct data sum;           //求和的结果

	int n;        
	int k;        //表示第k个数
    	
    //制造一个长度为N的序列 第k个数的值为k
	for(k=0;k<N;k++)
	{
		Squence[k]=k;
	}
	
	//进行DFT计算，并输出结果
	for(k=0;k<N;k++)
	{
		//对sum赋初值
		sum.real=0;               
	    sum.imag=0;

		//计算第k个点处的幅值
		for(n=0;n<N;n++)
		{
			one[n].real = cos(2*PI/N*k*n)*Squence[n];  //实部
			one[n].imag = -sin(2*PI/N*k*n)*Squence[n]; //虚部

			sum.real += one[n].real;                   //实部求和
			sum.imag += one[n].imag;                   //虚部求和
		}
		Amplitude[k] = sqrt(sum.real*sum.real+sum.imag*sum.imag);  //计算幅值 sqrt（a^2+b^2）
	  
		printf("%d %f\n",k,Amplitude[k]); //输出计算结果
	}
	return 0;
}

作图命令：

cd path                //打开文件路径
tcc dft.c -o dft.exe   //生成exe文件
dft.exe > dft.dat      //生成数据文件
gnuplot                //调用gnuplot
plot "dft.dat" wxt     //生成散点图

作图结果：