C语言实现DFT算法

一、什么是DFT？

在计算机机上实现信号的频谱分析及其他方面的处理工作时，对型号的要求是：在时域和频域都应是离散的，且都应是有限长。由于 $e(\pm j \frac{2\pi}{N} nk)$ 相对n和k都是以N为周期的，X(k)和x(n)为傅里叶变化对，即DFT。
DFT并不是一个新的傅里叶变换形式，它实际上来自于DFS，只不过仅在时域频域各取一个周期而已，由这一个周期作延拓。

二、DFT的作用

可以将信号从时域变换到频域，而且时域和频域都是离散的，通俗的说，可以求出一个信号由哪些正弦波叠加而成，求出的结果就是这些正弦波的幅度和相位。

三、DFT的实现

DFT公式：
$X(k)=\sum\limits^{N-1}_{n=0}x(n)e^{(-j \frac{2\pi}{N} kn)}$

X(k) — DFT变换后的结果

x(n) — 采样信号（实际中x(n)是实信号，即虚部为0）

1. 利用欧拉公式： $e^{j \theta}=cos\theta+jsin\theta$ 进行代换

$X(k)=\sum\limits^{N-1}_{n=0}x(n)cos(\frac{2\pi}{N} kn) -jx(n)sin(\frac{2\pi}{N} kn)$

实部： $real(k)=\sum\limits^{N-1}_{n=0}x(n)cos(\frac{2\pi}{N} kn)$
虚部： $imag(k)=\sum\limits^{N-1}_{n=0}-jx(n)sin(\frac{2\pi}{N} kn)$

2. 将sin和cos合并

$asin(x)+bcos(x)=\sqrt{a^2+b^2}sin(x+arctanb/a)$

最终利用计算机编程出来的算法就是计算出幅值

3. 源程序【DFT.c】

#include <stdio.h>
#include <stdlib.h>
#include <math.h>			 
#define PI	3.1415926535	

int N ;    //定于序列长度变量
double Input_Squence[100];	//输入的原始数据序列 
double Amplitude[100] ;   //存储幅值计算结果

typedef struct				//复数结构体,用于实现傅里叶运算 
{
	double real,imag;   
}complex;
complex Result_Point[100];		

void DFT_Calculate_Point(int k)
{
	int n = 0;
	complex Sum_Point;
	complex One_Point[N];
	Sum_Point.real = 0;
	Sum_Point.imag = 0;
	for(n=0; n<N; n++)
	{
		One_Point[n].real = cos(2*PI/N*k*n)*Input_Squence[n];  //运用欧拉公式把复数拆分成实部和虚部 
		One_Point[n].imag = -sin(2*PI/N*k*n)*Input_Squence[n];
		
		Sum_Point.real += One_Point[n].real;//对每一个点的实部求和
		Sum_Point.imag += One_Point[n].imag;
		//对每一个点的虚部求和
	}
	Result_Point[k].real = Sum_Point.real;
	Result_Point[k].imag = Sum_Point.imag;
}

void DFT_Calculate()
{
	int i = 0;
	for(i=0; i<N; i++)
	{
		DFT_Calculate_Point(i);
		Amplitude[i] = sqrt(Result_Point[i].real * Result_Point[i].real + Result_Point[i].imag * Result_Point[i].imag);  //计算幅值
	}
}

int main(int argc, char *argv[])
{
	N = atoi(argv[1]);  //命令行参数argv[1]为输入的幅值
	int i = 0;
	for(i=0; i<N; i++)//制造输入序列 
	{
		Input_Squence[i] = i;
	}
	DFT_Calculate(); //进行DFT计算 
	for(i=0; i<N; i++)
	{
		printf("%d\t%lf\n",i,Amplitude[i]); //输出计算结果
	}
	return 0;
}

四、计算结果验证

1. 运行【DFT.c】文件并通过gnuplot画图

Input_Squence={0，1，2，3，4，5，6，7，8，9，10，11，12，13，14，15}
N=16

2. 利用matlab计算DFT并画图

DFT函数定义：【DFT.m】

function [ Xk ] = DFT( xn,N ) %定于DFT算法
xxn=zeros(1,N);
xxn([1:length(xn)])=xn; %构造长度为N的序列
Xk=ones(1,N);
Wn=exp(-j*2*pi/N);
wn=ones(N,N);
for k=0:1:N-1   %DFT公式
    for n=0:1:N-1
        wn(k+1,n+1)=Wn^(k*n);
    end
end
Xk=(wn*xxn')';
end

利用DFT进行计算：

xn=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]; %输入的采样序列
N=16;  %序列长度
Xk3=DFT(xn,N);  %调用DFT
L=0:1:N-1;  %横轴长度
subplot(1,1,1);  %画图
stem(L,abs(Xk3));  %数据源
xlabel('k');  %设置横轴名称
ylabel('X(k)');   %设置纵轴名称
title('DFT N=16');

在这里插入图片描述