目錄

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

FFT修改

FFT-dit的算法討論

級(jí)的概念

$M=lo{g}_{2}N$,所以N點(diǎn)DFT可分成M級(jí)。

下圖為8點(diǎn)FFT時(shí)間抽取算法信號(hào)流程圖，方框分別為m=0,m=1,m=2;

碼位倒置

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

雷德算法： 對(duì)于自然順序（二進(jìn)制）我們是在低位加 1 得到下一位數(shù)，對(duì)于倒位序我們是在高位加 1 向低位進(jìn)位。比如已知一個(gè)倒位序數(shù)是J求其下一個(gè)倒位序數(shù)，N位總數(shù) ，把J與N/2比較若J<N/2則J的最高位為 0 ，把最高位置 1 ，就得到了J的下一個(gè)倒位序數(shù)；若J>=N/2則說明J的最高位為1 ，把最高位置0 ，比較次高位，若次高位為0 ，則把次高位置1，得到J的下一個(gè)倒位序，若次高位為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; }

結(jié)果：

$W^r$因子的分布

規(guī)律：第m級(jí)，$W_{2^{m+1}}^r$

在fft中，乘法主要來自旋轉(zhuǎn)因子。因?yàn)?span id="ozvdkddzhkzd" class="katex--inline">$W^r=cos(\frac{2\pi r}{N})?jsin(\frac{2\pi r}{N}),所以在對(duì)W^r$相乘時(shí)，必須產(chǎn)生相應(yīng)的正，余弦函數(shù)。編程時(shí)，產(chǎn)生正，余弦函數(shù)兩個(gè)方法，一個(gè)是在每一步直接產(chǎn)生，二是在程序開始前預(yù)先計(jì)算W^r存于一個(gè)數(shù)組中，等效一個(gè)正，余弦函數(shù)的“表”。

蝶形單元

第m級(jí)，有

$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)$

"組"的概念

每一級(jí)的N/2個(gè)蝶形單元可以分成若干組，每一組有著相同的結(jié)構(gòu)及W^r因子分布。

第m級(jí)的組數(shù)是$N/2^{m+1}$.

程序

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

B.FFT蝶形算法使用三重循環(huán)，下一層的數(shù)據(jù)都是由上一層計(jì)算得到的。

#include <stdio.h> #include<stdlib.h> #define pi 3.1415 #define N 64 typedef struct {double real,imag; } complex; //復(fù)數(shù)結(jié)構(gòu) complex fft_outa[100];//a的全部數(shù)據(jù) complex fft_outao[100];//a的單個(gè)數(shù)據(jù) complex fft_outb[100];//b的全部數(shù)據(jù) complex fft_outbo[100];//b的單個(gè)數(shù)據(jù) 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])//三層循環(huá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++) //第一層循環(huán){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);//旋轉(zhuǎn)因素變化速度//初始化旋轉(zhuǎn)因子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;//蝶形運(yùn)算的兩個(gè)輸入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++)//輸入信號(hào){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結(jié)果：

matlab的程序結(jié)果：

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

對(duì)比結(jié)果還是有差異。

FFT過程中，對(duì)C理解還不夠，對(duì)蝶形運(yùn)算算法還不理解，對(duì)三重循環(huán)還不熟悉 。

修改了，matlab的程序

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

得出和C一樣的結(jié)果

總結(jié)

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

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