設計需求

• 根據離散傅里葉變換的原始公式和自己編寫復數計算函數進行離散傅里葉變換
• 對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-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]);

總結

以上是生活随笔為你收集整理的C语言进行离散傅里叶DFT变换~MATLAB验证的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。