关于MATLAB实现以及C语言实现离散傅里叶变换（DFT）——傅里叶变换的学习与认识（二）

3.离散傅里叶变换（DFT）

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

用C语言实现DFT

设有限长序列的长度为N，对应的离散傅里叶变换定义为：
$f(n)=\sum_{k=0}^{N-1} f(k) e^{-j\frac{2\pi}{N}nk}$
我们用C语言实现DFT计算，实际就是求f(n)的值，也就是求信号的幅值。

根据欧拉公式
$e^{j\theta}=cos\theta + jsin\theta$
在c语言中复数由实部和虚部组成，代入DFT的计算式可以得到：
$f(n)=\sum_{k=0}^{N-1} f(k)cos(\frac{2\pi}{N}nk)-jsin(\frac{2\pi}{N}nk)$
实现代码

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

#define PI 3.14159

struct data
{
	double real,imag;
};

int main()
{
	int N=32;                  //序列长度
double Squence[100];       //序列
double Amplitude[100];     //幅值

struct data one[32];       //单个数据
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 \t%f\n",k,Amplitude[k]); //输出计算结果
}
return 0;
}

调用tcc与gnuplot编译与做图

在DFT.c的路径下

tcc DFT.c  //将源文件编译生成为可执行文件DFT.exe

DFT.exe > DFT.txt //通过重定向的方法，将执行结果产生的数据重新放到一个新的文件中

gnuplot //调用画图工具

plot[0:32] [0:120]  "DFT.txt" u 1:2 w p //作出该信号的离散傅里叶变换的点状图

运行结果

MATLAB实现DFT

实现代码

%计算长度为32点的DFT
N=32;%信号的长度
x=[1:N];%定义存放信号的数组
a=zeros(1,N);b=zeros(1,N);%构造长度为N的序列
for k=0:N-1
    for i=0:N-1
        a(k+1)=a(k+1)+x(i+1)*cos(2*pi*k*i/N);%DFT公式
        b(k+1)=b(k+1)+x(i+1)*sin(2*pi*k*i/N);
    end
c(k+1)=sqrt(a(k+1)^2+b(k+1)^2);%算出幅值
end

f=(0:1:N-1);
plot(f,c,'r'); %画图
axis([0 32 0 120]);
title('DFT');xlabel('k');ylabel('x(k)');

运行结果

结果对比

发现整体误差较小