設(shè)計(jì)需求

• 根據(jù)離散傅里葉變換的原始公式和自己編寫復(fù)數(shù)計(jì)算函數(shù)進(jìn)行離散傅里葉變換
• 對(duì)10000個(gè)點(diǎn)的加有噪聲或干凈的正弦波的數(shù)據(jù)進(jìn)行離散傅里葉變換，生成10000個(gè)點(diǎn)的復(fù)數(shù)數(shù)據(jù)序列到文本文件中。
• 數(shù)據(jù)格式為實(shí)部+虛部，用空格或逗號(hào)隔開。

實(shí)現(xiàn)思路

離散傅里葉變換的公式如下：
$$\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)$$
幅值計(jì)算：
$$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;//實(shí)部double imag;//虛部 }complex; //復(fù)數(shù)加法 void complex_add(complex a, complex b, complex *c) {c -> real = a.real + b.real;c -> imag = a.imag + b.imag; } //復(fù)數(shù)乘法 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) //定義旋轉(zhuǎn)因子 {//用歐拉公式分成實(shí)部和虛部 Wn->real=cos(2*pi*i*k/n1);Wn->imag=-sin(2*pi*i*k/n1); } int main(int argc, char* argv[]) {//檢查命令行參數(shù)正確與否if (argc != 4){printf("用法:程序名 輸入文件 輸出文件 采樣點(diǎn)數(shù)\n");printf("程序名：編譯后后的.exe名稱\n");printf("輸入文件：所需進(jìn)行計(jì)算的信號(hào)文件\n");printf("輸出文件：計(jì)算后輸出的文件名\n");printf("采樣點(diǎn)數(shù)：所需計(jì)算的采樣點(diǎn)數(shù)\n");return -1;}FILE* fp1 = fopen(argv[1], "r");int N = atoi(argv[2]);//定義采樣點(diǎn)數(shù)FILE* fp2 = fopen(argv[2], "w");int len = 0;double data = 0;//檢查文件能否被打開if (fp1 == NULL){printf("file error！\n");return -1;}//檢測(cè)文件行數(shù)len = data_len(fp1);complex x[len];//原序列complex X[len];//DFT后的序列double y[len];//幅值printf("data len:%d\n",len);rewind(fp1);//拷貝文件內(nèi)容char* str = (char*)malloc(sizeof(char) * len);for(int i = 0;i < len; i++)//把數(shù)據(jù)取出來{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++)//{//旋轉(zhuǎn)因子計(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++)//幅值計(jì)算// {// 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中編譯文件并運(yùn)行

打開result1.txt,查看數(shù)據(jù)

將數(shù)據(jù)導(dǎo)入gnuplot中繪制頻譜圖

在MATLAB中驗(yàn)證結(jié)果

MATLAB中的代碼如下

clear; close all; x = load('wave.txt');%讀取數(shù)據(jù) N = length(x);%檢查矩陣的長(zhǎng)度 a = zeros(N,1);%建立N行0矩陣，存放實(shí)部數(shù)據(jù) b = zeros(N,1);%建立N行0矩陣，存放虛部數(shù)據(jù) 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);%實(shí)部endc(k+1)=sqrt(a(k+1)^2+b(k+1)^2);%求幅度 end t=(0:1:N-1);%時(shí)間數(shù)據(jù) plot(t,c);%畫頻譜圖 axis([0,300,0,5000]);

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。