C语言进行离散傅里叶DFT变换~MATLAB验证

设计需求

根据离散傅里叶变换的原始公式和自己编写复数计算函数进行离散傅里叶变换
对10000个点的加有噪声或干净的正弦波的数据进行离散傅里叶变换，生成10000个点的复数数据序列到文本文件中。
数据格式为实部+虚部，用空格或逗号隔开。

实现思路

离散傅里叶变换的公式如下：
X(k)=∑n=0N−1x(n)exp⁡(−j2πNnk)=∑n=0N−1x(n)WNnk\begin{aligned} X(k)&=\sum_{n=0}^{N-1}x(n)\exp(-j\frac{2\pi}{N}nk)\\ &=\sum_{n=0}^{N-1}x(n)W_N^nk \end{aligned} X(k)=n=0∑N−1x(n)exp(−jN2πnk)=n=0∑N−1x(n)WNnk
带入欧拉公式
eix=cos⁡x+i(sin⁡x)e^{ix}=\cos x+i(\sin x)eix=cosx+i(sinx)
得到：
X(k)=∑n=0N−1x(n)cos⁡(2πNkn)−jx(n)sin⁡(2πNkn)X(k)=\sum_{n=0}^{N-1}x(n)\cos (\frac{2\pi}{N}kn)-jx(n)\sin (\frac{2\pi}{N}kn)X(k)=n=0∑N−1x(n)cos(N2πkn)−jx(n)sin(N2πkn)
幅值计算：
y(k)=a2+b2y(k)=\sqrt{a^2+b^2}y(k)=a2+b2

DFT源代码

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define pi 3.14159
typedef struct
{double real;//实部double imag;//虚部
}complex;
//复数加法
void complex_add(complex a, complex b, complex *c)
{c -> real = a.real + b.real;c -> imag = a.imag + b.imag;
}
//复数乘法
void complex_mul(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;
}
//求模
double complex_mod(complex a)
{double c = 0;c = pow(a.real * a.real + a.imag + a.imag, 0.5);return c;
}
int data_len(FILE*fp){int len = 0,c;while ((c = fgetc(fp)) != EOF)if (c == '\n'){len++;}return len;
}
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);
}
int main(int argc, char* argv[])
{//检查命令行参数正确与否if (argc != 4){printf("用法:程序名 输入文件 输出文件 采样点数\n");printf("程序名：编译后后的.exe名称\n");printf("输入文件：所需进行计算的信号文件\n");printf("输出文件：计算后输出的文件名\n");printf("采样点数：所需计算的采样点数\n");return -1;}FILE* fp1 = fopen(argv[1], "r");int N = atoi(argv[2]);//定义采样点数FILE* fp2 = fopen(argv[2], "w");int len = 0;double data = 0;//检查文件能否被打开if (fp1 == NULL){printf("file error！\n");return -1;}//检测文件行数len = data_len(fp1);complex x[len];//原序列complex X[len];//DFT后的序列double y[len];//幅值printf("data len:%d\n",len);rewind(fp1);//拷贝文件内容char* str = (char*)malloc(sizeof(char) * len);for(int i = 0;i < len; i++)//把数据取出来{fscanf(fp1, "%s", str);if (feof(fp1)) break;data = atof(str);x[i].real = data;x[i].imag = 0;}for(int j =0;j <len; j++)//{//旋转因子计算complex sum = {0,0};for(int n = 0;n < len; n++){complex wn,t;Wn(len,n,j,&wn);complex_mul(x[n],wn,&t);complex_add(sum,t,&sum);}X[j].real = sum.real;X[j].imag = sum.imag;fprintf(fp2,"%lf %lf\n",X[j].real,X[j].imag);}// for(int k=0;k < len; k++)//幅值计算// {// double t = k/len;// y[k] = sqrt(X[k].real * X[k].real// + X[k].imag * X[k].imag);// fprintf(fp2,"%lf %lf\n",t,y[k]);// }printf("success");fclose(fp1);fclose(fp2);return 0;
}

在tcc中编译文件并运行

打开result1.txt,查看数据

将数据导入gnuplot中绘制频谱图

在MATLAB中验证结果

MATLAB中的代码如下

clear;
close all;
x = load('wave.txt');%读取数据
N = length(x);%检查矩阵的长度
a = zeros(N,1);%建立N行0矩阵，存放实部数据
b = zeros(N,1);%建立N行0矩阵，存放虚部数据
for k=0:N-1for i=0:N-1b(k+1)=b(k+1)+x(i+1)*sin((2*pi*k*i)/N);%虚部a(k+1)=a(k+1)+x(i+1)*cos((2*pi*k*i)/N);%实部endc(k+1)=sqrt(a(k+1)^2+b(k+1)^2);%求幅度
end
t=(0:1:N-1);%时间数据
plot(t,c);%画频谱图
axis([0,300,0,5000]);