目錄

• FFT修改
• • FFT-dit的算法討論
  • 程序
  • • c結果：
    • matlab的程序結果：
    • 修改了，matlab的程序

FFT修改

FFT-dit的算法討論

級的概念

$M=lo{g}_{2}N$,所以N點DFT可分成M級。

下圖為8點FFT時間抽取算法信號流程圖，方框分別為m=0,m=1,m=2;

碼位倒置

交換后輸出序列x(k)依照正序排列，但輸入序列x(n)的次序不再是原來的自然次序，這是由于將x(n)按奇,偶分開產生的。

雷德算法： 對于自然順序（二進制）我們是在低位加 1 得到下一位數，對于倒位序我們是在高位加 1 向低位進位。比如已知一個倒位序數是J求其下一個倒位序數，N位總數 ，把J與N/2比較若J<N/2則J的最高位為 0 ，把最高位置 1 ，就得到了J的下一個倒位序數；若J>=N/2則說明J的最高位為1 ，把最高位置0 ，比較次高位，若次高位為0 ，則把次高位置1，得到J的下一個倒位序，若次高位為1 ， 則把次高位置0,以此類推…

參考：https://blog.csdn.net/corcplusplusorjava/article/details/17119567

Rader參考程序：

#include <stdio.h> #include <stdlib.h> #include <math.h> int main(void) {int array[8]={0,1,2,3,4,5,6,7};int i,j,k;int N = 8;int temp;j = 0;for(i = 0; i < N -1; i ++){if(i < j){temp = array[i];array[i] = array[j];array[j] = temp;}k = N >> 1;while( k <= j){j = j - k;k >>= 1;}j = j + k;}for( i = 0; i < N; i ++)printf("%d %d\n",i,array[i]);printf("\n");return 0; }

結果：

$W^r$因子的分布

規律：第m級，$W_{2^{m+1}}^r$

在fft中，乘法主要來自旋轉因子。因為$W^r=cos(\frac{2\pi r}{N})?jsin(\frac{2\pi r}{N}),所以在對W^r$相乘時，必須產生相應的正，余弦函數。編程時，產生正，余弦函數兩個方法，一個是在每一步直接產生，二是在程序開始前預先計算W^r存于一個數組中，等效一個正，余弦函數的“表”。

蝶形單元

第m級，有

$X_{m+1}(p)=x_m(p)+W_N^r{X}_m(q)$

$X_{m+1}(q)=x_m(p)?W_N^r{X}_m(q)$

"組"的概念

每一級的N/2個蝶形單元可以分成若干組，每一組有著相同的結構及W^r因子分布。

第m級的組數是$N/2^{m+1}$.

程序

A.1<<i是把1左移i位，每次左移以為就是乘以2，所以1<<i的結果是1乘以2的i次方 i<<1就是把i左移一位，即i乘以2，假如i=5,最后結果就是5*2=10

B.FFT蝶形算法使用三重循環，下一層的數據都是由上一層計算得到的。

#include <stdio.h> #include<stdlib.h> #define pi 3.1415 #define N 64 typedef struct {double real,imag; } complex; //復數結構 complex fft_outa[100];//a的全部數據 complex fft_outao[100];//a的單個數據 complex fft_outb[100];//b的全部數據 complex fft_outbo[100];//b的單個數據 double x_r[N],x_i[N]; double y_r[N],y_[N];void Rader(double x_r[N],double x_i[N])//雷德算法排序 {int i,j,k;j=0;double t_r,t_i;for(i=0; i<N-1; i++){if(i<j){t_r=x_r[i];t_i=x_i[i];x_r[i]=x_r[j];x_i[i]=x_i[j];x_r[j]=t_r;x_i[j]=t_i;}k=N>>1;while(k<=j){j=j-k;k>>=1;}j=j+k;} } void dit(double x_r[N],double x_i[N])//三層循環 {int m,r,p,q;int step,factor_step;int a,b,i;double t_real,t_imag;double w_re,w_im;for(m=0; m<log2(N); m++) //第一層循環{a=1;for(i=0;i<m+1;i++)a=a*2;b=a/2;w_re=cos(pi/b);w_im=-sin(pi/b);step=1<<(m+1);// factor_step=N>>(m+1);//旋轉因素變化速度//初始化旋轉因子double factor_real=1.0;double factor_imag=0.0;for(r=0; r<step/2; r++){// w_re=cos(2*pi*(r+1)*factor_step/N);// w_im=-sin(2*pi*(r+1)*factor_step/N);for(p=r; p<N; p+=step){q=p+step/2;//蝶形運算的兩個輸入t_real=x_r[q]*factor_real-x_i[q]*factor_imag;t_imag=x_r[q]*factor_imag+x_i[q]*factor_real;x_r[q]=x_r[p]-t_real;x_i[q]=x_i[p]-t_imag;x_r[p]=x_r[p]+t_real;x_i[p]=x_i[p]+t_imag;}t_real=factor_real*w_re-factor_imag*w_im;t_imag=factor_real*w_im+factor_imag*w_re;factor_real=t_real;factor_imag=t_imag;}}}int main() {int n,i;double y[N];for(n=0; n<N; n++)//輸入信號{x_r[n]=cos(n*pi/6);x_i[n]=0;}Rader(x_r,x_i);dit(x_r,x_i);for(i=0; i<N; i++)//輸出fft{y[i]=sqrt(x_r[i]*x_r[i]+x_i[i]*x_i[i]);printf("%d %f\n",i,y[i]);} }

c結果：

matlab的程序結果：

N=64; n=0:1:N; x=cos(n*pi/6); X=fft(x); figure stem(n,abs(X))

對比結果還是有差異。

FFT過程中，對C理解還不夠，對蝶形運算算法還不理解，對三重循環還不熟悉 。

修改了，matlab的程序

>> n=0:1:63; >> x=cos(n*pi/6); >> y=fft(x); >> figure >> stem(n,abs(y))

得出和C一樣的結果

總結

以上是生活随笔為你收集整理的FFT C语言 修改了matlab的全部內容，希望文章能夠幫你解決所遇到的問題。

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