fft

来源：互联网发布：mp4下载视频软件编辑：程序博客网时间：2024/06/06 12:35

#include <stdio.h> #include <stdlib.h> #include <math.h> //const float PI = 3.1416; inline void swap (float &a, float &b) { float t; t = a; a = b; b = t; } void bitrp (float xreal [], float ximag [], int n) { //位反转置换Bit-reversal Permutation int i, j, a, b, p; for (i = 1, p = 0; i < n; i *= 2) { p ++;} for (i = 0; i < n; i ++) { a = i; b = 0; for (j = 0; j < p; j ++) { b = (b << 1) + (a & 1); // b = b * 2 + a % 2; a >>= 1; // a = a / 2; } if ( b > i) { swap (xreal [i], xreal [b]); swap (ximag [i], ximag [b]); } } }void FFT(float xreal[], float ximag[], int n) { //快速傅立叶变换，将复数x 变换后仍保存在x 中，xreal, ximag 分别是x 的实部和虚部float wreal[128];float wimag[128];float treal,timag,ureal,uimag,arg;int m, k, j, t, index1, index2; bitrp (xreal, ximag, n); // 计算1 的前n / 2 个n 次方根的共轭复数W'j = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1 arg = - 2 * PI / n; treal = cos (arg); timag = sin (arg); wreal [0] = 1.0; wimag [0] = 0.0; for (j = 1; j < n / 2; j ++) { wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag; wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal; } for (m = 2; m <= n; m *= 2)  { for (k = 0; k < n; k += m){ for (j = 0; j < m / 2; j ++) { index1 = k + j;index2 = index1 + m / 2; t = n * j / m; // 旋转因子w 的实部在wreal [] 中的下标为t treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2]; timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2]; ureal = xreal [index1]; uimag = ximag [index1]; xreal [index1] = ureal + treal; ximag [index1] = uimag + timag; xreal [index2] = ureal - treal; ximag [index2] = uimag - timag; } } } } // // void IFFT (float xreal [], float ximag [], int n) // { // // 快速傅立叶逆变换// float wreal [128], wimag [128], treal, timag, ureal, uimag, arg; // int m, k, j, t, index1, index2; // bitrp (xreal, ximag, n); // // 计算1 的前n / 2 个n 次方根Wj = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1 // // arg = 2 * PI / n; // treal = cos (arg); // timag = sin (arg); // wreal [0] = 1.0; // wimag [0] = 0.0; // for (j = 1; j < n / 2; j ++) // { // wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag; // wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal; // } // for (m = 2; m <= n; m *= 2) // { // for (k = 0; k < n; k += m)  // { // for (j = 0; j < m / 2; j ++)// { // index1 = k + j; // index2 = index1 + m / 2; // t = n * j / m; // // 旋转因子w 的实部在wreal [] 中的下标为t // treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];// timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2]; // ureal = xreal [index1]; // uimag = ximag [index1];// xreal [index1] = ureal + treal;// ximag [index1] = uimag + timag; // xreal [index2] = ureal - treal; // ximag [index2] = uimag - timag; // } // } // } // for (j=0; j < n; j ++) // { // xreal [j] /= n; // ximag [j] /= n; // } // }

    #include "math.h"      #include "fft.h"      //精度0.0001弧度            void conjugate_complex(int n,complex in[],complex out[])      {        int i = 0;        for(i=0;i<n;i++)        {          out[i].imag = -in[i].imag;          out[i].real = in[i].real;        }       }            void c_abs(complex f[],float out[],int n)      {        int i = 0;        float t;        for(i=0;i<n;i++)        {          t = f[i].real * f[i].real + f[i].imag * f[i].imag;          out[i] = sqrt(t);        }       }                  void c_plus(complex a,complex b,complex *c)      {        c->real = a.real + b.real;        c->imag = a.imag + b.imag;      }            void c_sub(complex a,complex b,complex *c)      {        c->real = a.real - b.real;        c->imag = a.imag - b.imag;       }            void c_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;         }            void c_div(complex a,complex b,complex *c)      {        c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);        c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);      }            #define SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr            void Wn_i(int n,int i,complex *Wn,char flag)      {        Wn->real = cos(2*PI*i/n);        if(flag == 1)        Wn->imag = -sin(2*PI*i/n);        else if(flag == 0)        Wn->imag = -sin(2*PI*i/n);      }            //傅里叶变化      void fft(int N,complex f[])      {        complex t,wn;//中间变量        int i,j,k,m,n,l,r,M;        int la,lb,lc;        /*----计算分解的级数M=log2(N)----*/        for(i=N,M=1;(i=i/2)!=1;M++);         /*----按照倒位序重新排列原信号----*/        for(i=1,j=N/2;i<=N-2;i++)        {          if(i<j)          {            t=f[j];            f[j]=f[i];            f[i]=t;          }          k=N/2;          while(k<=j)          {            j=j-k;            k=k/2;          }          j=j+k;        }              /*----FFT算法----*/        for(m=1;m<=M;m++)        {          la=pow(2,m); //la=2^m代表第m级每个分组所含节点数               lb=la/2;    //lb代表第m级每个分组所含碟形单元数                       //同时它也表示每个碟形单元上下节点之间的距离          /*----碟形运算----*/          for(l=1;l<=lb;l++)          {            r=(l-1)*pow(2,M-m);               for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la            {              lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号                   Wn_i(N,r,&wn,1);//wn=Wnr              c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算              c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr              c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr            }          }        }      }            //傅里叶逆变换      void ifft(int N,complex f[])      {        int i=0;        conjugate_complex(N,f,f);        fft(N,f);        conjugate_complex(N,f,f);        for(i=0;i<N;i++)        {          f[i].imag = (f[i].imag)/N;          f[i].real = (f[i].real)/N;        }      }

matlab仿真

%例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。clf;fs=100;N=128;   %采样频率和数据点数n=0:N-1;t=n/fs;   %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N);    %对信号进行快速Fourier变换mag=abs(y);     %求得Fourier变换后的振幅f=n*fs/N;    %频率序列subplot(2,2,1),plot(f,mag);   %绘出随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;%对信号采样数据为1024点的处理fs=100;N=1024;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N);   %对信号进行快速Fourier变换mag=abs(y);   %求取Fourier变换的振幅f=n*fs/N;subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=1024');grid on;subplot(2,2,4)plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=1024');grid on;

DFT

#include <stdio.h>  #include <math.h> #include "fft.h" #define N   SAMPLEPOINTS              //序列长度        typedef float ElementType; //原始数据序列的数据类型,可以在这里设置   typedef struct              //复数结构体,用于实现傅里叶运算   {      ElementType real,imag;  }complex;  complex dataResult[N];      //傅里叶运算的频域结果序列的值(复数)  ElementType dataSource[N];  //输入的原始数据序列   ElementType dataFinualResult[N]; //最终频域序列结果(dataResult[N]的模,实数)  void FFT_Calculate_OneNode(int k)//计算频域上一个点的DFT值   {      int n = 0;      complex ResultThisNode;      complex part[N];      ResultThisNode.real = 0;      ResultThisNode.imag = 0;      for(n=0; n<N; n++)      {          part[n].real = cos(2*PI/N*k*n)*dataSource[n];//运用欧拉公式把复数拆分成实部和虚部           part[n].imag = sin(2*PI/N*k*n)*dataSource[n];          ResultThisNode.real += part[n].real;          ResultThisNode.imag += part[n].imag;      }      dataResult[k].real = ResultThisNode.real;      dataResult[k].imag = ResultThisNode.imag;  }  void FFT_Calculate()//对输入序列全部计算DFT  {      int i = 0;      for(i=0; i<N; i++)      {          FFT_Calculate_OneNode(i);          dataFinualResult[i] = sqrt(dataResult[i].real * dataResult[i].real + dataResult[i].imag * dataResult[i].imag);      }  }  int testdftmain(float * input,float * output)  {      int i = 0;      for(i=0; i<N; i++)//制造输入序列       {          dataSource[i] = input[i];              }      FFT_Calculate();//进行DFT计算           for(i=0; i<N; i++)  {output[i]=dataFinualResult[i];}      return 0;  }

0 0