3.离散傅里叶变换(DFT)
离散傅里叶变换指的是傅里叶变换在时域和频域上都呈离散的形式。在形式上,时域和频域上的序列都是有限长的,而实际上这两组序列都应当被认为是离散周期信号的主值序列。对于有限长的离散信号作DFT,可以把它看作是离散周期信号的一个周期,应进行周期延拓之后再进行计算。
用C语言实现DFT
设有限长序列的长度为N,对应的离散傅里叶变换定义为:
我们用C语言实现DFT计算,实际就是求f(n)的值,也就是求信号的幅值。
根据欧拉公式
在c语言中复数由实部和虚部组成,代入DFT的计算式可以得到:
实现代码
#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)');
运行结果
结果对比
发现整体误差较小
转载:https://blog.csdn.net/qq_40644291/article/details/102490851
查看评论