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

设计需求

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

实现思路

离散傅里叶变换的公式如下：
$$\begin{array}{rl}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^{n}k\end{array}$$
带入欧拉公式
$$e^{ix}=\cos x+i(\sin x)$$
得到：
$$X(k)=\sum _{n=0}^{N-1}x(n)\cos (\frac{2\pi }{N}kn)-jx(n)\sin (\frac{2\pi }{N}kn)$$
幅值计算：
$$y(k)=\sqrt{a^2+b^2}$$

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-1
    for i=0:N-1
        b(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);%实部
    end
    c(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]);