徐士良老師編寫的c語言算法程序下載鏈接：https://pan.baidu.com/s/1zDV6iLeYeXmZaoZlP4yRAA 提取碼：8opo

一：徐士良老師編寫的FFT程序比較
1.FFT源程序【kfft.c】
抽樣點為：64
CLOCKS_PER_SEC 單位1000ms

#include "math.h"void kfft(pr,pi,n,k,fr,fi)int n,k;double pr[],pi[],fr[],fi[];{ int it,m,is,i,j,nv,l0;double p,q,s,vr,vi,poddr,poddi;for (it=0; it<=n-1; it++) //將pr[0]和pi[0]循環賦值給fr[]和fi[]{ m=it; is=0;for(i=0; i<=k-1; i++){ j=m/2; is=2*is+(m-2*j); m=j;}fr[it]=pr[is]; fi[it]=pi[is];}pr[0]=1.0; pi[0]=0.0;p=6.283185306/(1.0*n);pr[1]=cos(p); //將w=e^-j2pi/n用歐拉公式表示pi[1]=-sin(p);for (i=2; i<=n-1; i++) //計算pr[]{ p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);pr[i]=p-q; pi[i]=s-p-q;}for (it=0; it<=n-2; it=it+2) { vr=fr[it]; vi=fi[it];fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];}m=n/2; nv=2;for (l0=k-2; l0>=0; l0--) //蝴蝶操作{ m=m/2; nv=2*nv;for (it=0; it<=(m-1)*nv; it=it+nv)for (j=0; j<=(nv/2)-1; j++){ p=pr[m*j]*fr[it+j+nv/2];q=pi[m*j]*fi[it+j+nv/2];s=pr[m*j]+pi[m*j];s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);poddr=p-q; poddi=s-p-q;fr[it+j+nv/2]=fr[it+j]-poddr;fi[it+j+nv/2]=fi[it+j]-poddi;fr[it+j]=fr[it+j]+poddr;fi[it+j]=fi[it+j]+poddi;}}for (i=0; i<=n-1; i++){ pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]); //幅值計算} return;}

2.源程序【FFT2.c】
輸入信號：0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <D:\temp\kfft.c> #include <time.h>#define PI 3.1415926 #define N 64typedef struct {double real;double imag; }complex;void Wn(int ,int ,int ,complex *);//定義旋轉因子 void c_jiafa(complex ,complex ,complex *); //復數加法運算 void c_chengfa(complex ,complex ,complex *); //復數乘法運算 //void kfft(preal,pimag,n,k,freal,fimag);int main() { complex f[N], A[N], a[N]; double y[N];int i,j;double pr[64],pi[64],fr[64],fi[64],t[64];clock_t begin, end;double cost1 ,cost2, persent;for (i=0; i<=63; i++) //生成輸入信號{ t[i] = i*0.001;pr[i]=0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i); pi[i]=0.0; //0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i)}begin = clock(); //開始記錄時間kfft(pr,pi,64,6,fr,fi); //調用FFT函數end = clock(); //結束記錄時間cost1 = (double)(end - begin) / CLOCKS_PER_SEC;printf("*****************************FFT變換輸出*****************************\n");for (i=0; i<64; i++){ printf("%d\t%lf\n",i,pr[i]); //輸出結果}printf("*****************************DFT變換輸出*****************************\n");for(int i=0;i<N;i++) //提供初始信號 {a[i].real=0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i);//0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)a[i].imag=0;}begin = clock(); //開始記錄時間for(int k=0;k<N;k++){//int n;complex sum={0,0};for(int n=0; n<N;n++){complex wn,t;Wn(N,n,k,&wn);c_chengfa(a[n],wn,&t);c_jiafa(sum,t,&sum);}A[k].real=sum.real;A[k].imag=sum.imag;} end = clock(); //結束記錄時間cost2 = (double)(end - begin) / CLOCKS_PER_SEC;for(int i=0; i<N;i++)//幅值計算 {y[i]=sqrt(A[i].real*A[i].real+A[i].imag*A[i].imag);printf("%d\t%lf\n",i,y[i]);}printf("在抽樣點為%d的情況下，FFT運算一共花費時間為：%lf秒\n",N, cost1); //CLOCKS_PER_SEC 單位1000msprintf("在抽樣點為%d的情況下，DFT運算一共花費時間為：%lf秒\n",N, cost2); //CLOCKS_PER_SEC 單位1000ms//persent = (double)(cost1 / cost2);//printf("在抽樣點為%d的情況下，FFT算法運算時間為DFT算法運算時間的%lf倍\n", N, persent);system("pause");return 0; }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); }void c_jiafa(complex a,complex b,complex *c) //復數加法運算 {c->real=a.real+b.real; c->imag=a.imag+b.imag; }void c_chengfa(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; }

3.TCC下編譯

4.結果顯示

二、自己理解得FFT程序比較(有不足之處還望為我指出，謝謝大家了)
1.程序源碼
輸入信號：0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)
抽樣點為：64

#include <stdio.h> #include <math.h> #include <time.h>#define N 64 //設置抽樣點數 #define PI 3.1514926 //定義圓周率typedef struct //定義復數結構體變量 {double real;double imag; }complex;void c_jiafa(complex, complex, complex *); //復數加運算 void c_jianfa(complex, complex, complex *); //復數減運算 void c_chengfa(complex, complex, complex *); //復數乘運算 void Wn_i(int, int, complex *); //FFT旋轉因子 void Wn_ik(int, int, int, complex *); //DFT旋轉因子int main() {complex f[N], A[N], a[N]; //f[N]為輸出的快速傅里葉變換序列 A[N]為DFT變換輸出序列 a[N]為初始序列int LH, K, J, B, L, k, N1, P, M, K1; //L表示第L級蝶形 p旋轉因子指數 B兩序列間隔點數k 第k個序列double T;double y[N];double pr[N];clock_t begin, end;double cost1, cost2, persent; //定義時間double t_r,t_i;M =6;LH = N / 2;J = LH;N1 = N - 1;for (int i = 0; i < N; i++) //為DFT運算提供初始序列{a[i].real =0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i) ;a[i].imag = 0;}printf("******************************* 級數為：%d *******************************\n", M);//.......................................................................................................................for (int I = 1; I < N1; I++) //定義倒序序列函數{if (I < J){t_r=a[I].real;t_i=a[I].imag;a[I].real=a[J].real;a[I].imag=a[J].imag;a[J].real=t_r;a[J].imag=t_i;}K = LH;if (J >= K){do{J = J - K;K = K / 2;} while (J >= K);}J = J + K;}for (int i = 0; i < N; i++) //將序列賦給結構體函數進行運算{f[i].real = 0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i);f[i].imag = 0;}begin = clock(); //開始記錄時間for (L = 1; L <= M; L++) //FFT運算{B = (int)(pow(2, L - 1));for (J = 0; J < B; J++){P = (int)(J*pow(2, M - L));for (k = J; k < N; k = (int)(k + pow(2, L))){K1 = k + B;complex wn, t;Wn_i(N, P, &wn);c_chengfa(f[K1], wn, &t); //。。。。。。。。。。。。c_jianfa(f[k], t, &(f[K1])); //蝶形運算c_jiafa(f[k], t, &(f[k])); //。。。。。。。。。。。。}}}end = clock(); //結束記錄時間cost1 = (double)(end - begin) / CLOCKS_PER_SEC;//.......................................................................................................................printf("*****************************快速傅里葉變換輸出*****************************\n");for (int i=0; i<=N-1; i++){ pr[i]=sqrt(f[i].real*f[i].real+f[i].imag*f[i].imag); //幅值計算printf("%d\t%lf\n", i, pr[i] );}//.......................................................................................................................begin = clock(); //開始記錄時間for (int k = 0; k < N; k++) //DFT運算{complex sum = { 0, 0 };for (int n = 0; n < N; n++){complex wn, t;Wn_ik(N, n, k, &wn);c_chengfa(a[n], wn, &t);c_jiafa(sum, t, &sum);}A[k].real = sum.real;A[k].imag = sum.imag;}end = clock(); //結束記錄時間cost2 = (double)(end - begin) / CLOCKS_PER_SEC;//.......................................................................................................................printf("*****************************DFT變換輸出*****************************\n");for (int i = 0; i < N; i++) //DFT變換輸出{y[i]=sqrt(A[i].real*A[i].real+A[i].imag*A[i].imag);printf("%d\t%lf\n",i,y[i]);//printf("%d\t%lf\n", i, A[i].real, A[i].imag);}//.......................................................................................................................printf("FFT算法一共花費時間為：%lf秒\n", cost1); //CLOCKS_PER_SEC 單位1000msprintf("DFT算法一共花費時間為：%lf秒\n", cost2); //CLOCKS_PER_SEC 單位1000mspersent = (float)(cost1 / cost2);printf("在抽樣點為%d的情況下，FFT算法運算時間為DFT算法運算時間的%f倍\n", N, persent);return 0; }void c_jiafa(complex a, complex b, complex *c) //復數加法 {c->real = a.real + b.real;c->imag = a.imag + b.imag; } void c_jianfa(complex a, complex b, complex *c) //復數減法 {c->real = a.real - b.real;c->imag = a.imag - b.imag; } void c_chengfa(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; }void Wn_i(int n1, int i, complex *Wn) //定義FFT旋轉因子 {Wn->real = cos(2 * PI*i / n1);Wn->imag = -sin(2 * PI*i / n1); }void Wn_ik(int n1, int i, int k, complex *Wn) //定義DFT旋轉因子 {Wn->real = cos(2 * PI*i*k / n1);Wn->imag = -sin(2 * PI*i*k / n1); }

2.TCC編譯

3.結果顯示

三、總結
FFT是DFT的快速算法，可以節省大量的計算時間，快速傅里葉變換（FFT）是一種能在O(nlogn)的時間內將一個多項式轉換成它的點值表示的算法

總結

以上是生活随笔為你收集整理的对比DFT程序与FFT程序的效率的全部內容，希望文章能夠幫你解決所遇到的問題。

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