STM32--任意点数FFT变换程序

1万 · 2015-12-29 22:27

#include "math.h"

#define FFT_N 256                //----N点FFT变换

#define FFT_M 8                //----N等于2的M次方



u16 weizhi[FFT_N]; 

u16 wei1[FFT_N]; 

float temdat1[FFT_N];

float temdat2[FFT_N];

float tem_shi[FFT_N];

float tem_xu[FFT_N];  

float out_shi[FFT_N];

float out_xu[FFT_N];  

float Pix=3.1415926; 

float P=0;

float v1,v2;

float v3,v4,v5,v6; 







//----------------------------------------------------------------------------------

//--- 函数        

//f=(double)i/100; //-- i的取值为 0-255，步进为1            

//ddt=sin(2*pi*20*f);        

//--正弦函数，频率20Hz//--以fs=100Hz的采样频率采样，采样N=256个点，进行256点FFT运算

//----在数据点位x=51时--频率量的值达到最大点此时的频率量为 f=(x*fs)/(N-1)=20Hz

//----------------------------------------------------------------------------------

float ddt[FFT_N]=

{



0.000000,0.951057,0.587785,-0.587785,-0.951057,-0.000000,0.951056,0.587785,

-0.587785,-0.951057,-0.000000,0.951056,0.587785,-0.587785,-0.951057,-0.000000,

0.951056,0.587786,-0.587785,-0.951057,-0.000000,0.951056,0.587786,-0.587785,

-0.951057,-0.000001,0.951056,0.587786,-0.587785,-0.951057,-0.000001,0.951056,

0.587786,-0.587785,-0.951057,-0.000001,0.951056,0.587786,-0.587785,-0.951057,

-0.000001,0.951056,0.587786,-0.587785,-0.951057,-0.000001,0.951056,0.587786,

-0.587784,-0.951057,-0.000001,0.951056,0.587786,-0.587784,-0.951057,-0.000001,

0.951056,0.587786,-0.587784,-0.951057,-0.000001,0.951056,0.587786,-0.587784,

-0.951057,-0.000001,0.951056,0.587786,-0.587784,-0.951057,-0.000002,0.951056,

0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587784,-0.951057,

-0.000002,0.951056,0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,

-0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587784,-0.951057,-0.000002,

0.951056,0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587783,

-0.951057,-0.000002,0.951056,0.587787,-0.587783,-0.951057,-0.000002,0.951056,

0.587787,-0.587783,-0.951057,-0.000002,0.951056,0.587787,-0.587783,-0.951057,

-0.000003,0.951056,0.587787,-0.587783,-0.951057,-0.000003,0.951056,0.587787,

-0.587783,-0.951057,-0.000003,0.951056,0.587788,-0.587783,-0.951057,-0.000003,

0.951056,0.587788,-0.587783,-0.951057,-0.000003,0.951056,0.587788,-0.587783,

-0.951057,-0.000003,0.951056,0.587788,-0.587783,-0.951058,-0.000003,0.951056,

0.587788,-0.587783,-0.951058,-0.000003,0.951055,0.587788,-0.587783,-0.951058,

-0.000003,0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,0.587788,

-0.587782,-0.951058,-0.000004,0.951055,0.587788,-0.587782,-0.951058,-0.000004,

0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,0.587788,-0.587782,

-0.951058,-0.000004,0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,

0.587789,-0.587782,-0.951058,-0.000004,0.951055,0.587789,-0.587782,-0.951058,

-0.000004,0.951055,0.587789,-0.587782,-0.951058,-0.000004,0.951055,0.587789,

-0.587782,-0.951058,-0.000005,0.951055,0.587789,-0.587782,-0.951058,-0.000005,

0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587789,-0.587781,

-0.951058,-0.000005,0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,

0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587789,-0.587781,-0.951058,

-0.000005,0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587790,

-0.587781,-0.951058,-0.000005,0.951055,0.587790,-0.587781,-0.951058,-0.000005,





};





//------------------------------------------------------------------------------

//-- 复数 求模  ---------------------------------------------------------------- 

//------------------------------------------------------------------------------

float DAT_sqrt(float A1,float A2)

{

   u8 i=0; 

   float x;

         float y;

         x=A1*A1;

   y=A2*A2;

         y=x+y;

   if(y<1) x=1;

   else  x=y/2;  

   for(i=0;i<15;i++)

   {    

     if(x*x<=y)

                                break;

     else

        x=((y/x)+x)/2; //牛顿迭代法  

   }

   return x;

}





//-----------------------------------------------------------------------------

//---快速傅里叶变换 N个点变换--M层 蝶形计算 

//----------------------------------------------------------------------------- 

void FFT(float in_dat[],float *jieguo)

{

   u16 L,B; 

   u16 i=0,j=0;

   u16 j1,j2=0;

         u16 N=FFT_N;

         u16 M=FFT_M;

        

//-----------------------------位置置换 b2b1b0---b0b1b2--做M次左移 

   for(i=0;i<N;i++)

   { 

     weizhi[i]=i;

     wei1[i]=0;

   }

   for(i=0;i<N;i++)

   {

      for(j=0;j<M;j++)

      {

        wei1[i]=wei1[i]<<1;

        if(weizhi[i]&(0x01<<j))

        wei1[i]|=0x01; 

        else

        wei1[i]&=0xfffe;              

      } 

      weizhi[i]=wei1[i];                      

   } 

   for(i=0;i<N;i++)

   temdat1[i]=in_dat[i];

   for(i=0;i<N;i++)

   temdat2[i]=temdat1[weizhi[i]];

//------------------------------------------------------------------------------   

//------------------------------------------------------------------------------

//----变换位置的数据存在 temdat2[]中 

//---- 欧拉公式 e^(-jx)=cos(x)-jsin(x);     

//------蝶形计算----------------------------------------------------------------

//---B表示蝶形中两点的距离---L表示第L级蝶形运算--j表示的j行的计算 

//------------------------------------------------------------------------------ 

  for(i=0;i<N;i++) 

  {

   out_shi[i]=temdat2[i];

   out_xu[i]=0.0;

  }

  B=1; 

                                

  for(L=1;L<=M;L++)

  {

   for(j=0;j<N;j++)

   {

     j1=1<<L;

     j2=1<<(L-1);

     if((j%j1)<j2) 

     { 

      P=j*N;

      for(i=0;i<(L-1);i++)

      P=P/2;

      P=P*Pix;

      P=P/N;

      v3=cos(P);

      v4=0-sin(P);

      v5=(out_shi[j+B]*v3);

      v6=(out_xu[j+B]*v4);   

      v1=v5-v6;        

      v5=(out_shi[j+B]*v4);

      v6=(out_xu[j+B]*v3);

      v2=v5+v6;                               

      tem_shi[j]  =out_shi[j]+v1;         //--计算实部                  

      tem_xu[j]   =out_xu[j]+v2;          //--计算虚部  

      tem_shi[j+B]=out_shi[j]-v1;         //--计算实部       

      tem_xu[j+B] =out_xu[j]-v2;          //--计算虚部  

      out_shi[j]=tem_shi[j];

      out_xu[j] =tem_xu[j];

      out_shi[j+B]=tem_shi[j+B];

      out_xu[j+B] =tem_xu[j+B];                                    

     }            

   }

   B=B*2;                 

  } 

//-------------------------------------------------------

//---FFT 变换后的复数取模--------------------------------

//-------------------------------------------------------

        for(i=0;i<N;i++) 

        {

                *jieguo=DAT_sqrt(out_shi[i],out_xu[i]);

                jieguo++;

        } 

//-------------------------------------------------------   

}







//----------------------------------------------

//--变换结果如下

//----------------------------------------------

void shujuchuli(void)

{



        FFT(ddt,ddt); 





}









MATLAB软件计算程序



%---------------------------------------------------------

clf;

fs=100;N=256;   %采样频率和数据点数



n=0:N-1;t=n/fs;   %时间序列



x=sin(2*pi*20*t); %信号



y=fft(x,N);    %对信号进行快速Fourier变换

mag=abs(y);     %求得Fourier变换后的振幅



f=n*fs/N;    %频率序列



subplot(2,2,3),plot(n,x);   %绘出随频率变化的振幅

xlabel('t/Hz');

ylabel('振幅');title('N=256');grid on;





subplot(2,2,1),plot(f,mag);   %绘出随频率变化的振幅

xlabel('频率/Hz');

ylabel('振幅');title('N=256');grid on;





subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅

xlabel('频率/Hz');

ylabel('振幅');title('N=256');grid on;



%---------------------------------------------------------

1万 · 2015-12-29 22:47

一，如果对信号进行同样点数N的FFT变换，采样频率fs越高，则可以分析越高频的信号；与此同时，采样频率越低，对于低频信号的频谱分辨率则越好。

二，假设采样点不在正弦信号的波峰、波谷、以及0电压处，频谱则会产生泄露（leakage）。

三，对于同样的采样率fs，提高FFT的点数N，则可提高频谱的分辨率。

四，如果采样频率fs小于2倍信号频率2*fs（奈圭斯特定理），则频谱分析结果会出错。

1万 · 2015-12-29 22:47

复制

%清除命令窗口及变量

clc;

clear all;



%输入f、N、T、是否补零（补几个零）

f=input('Input frequency of the signal: f\n');

N=input('Input number of pointsl: N\n');

T=input('Input sampling time: T\n');

flag=input('Add zero too sampling signal or not? yes=1 no=0\n');

if(flag)

    ZeroNum=input('Input nmber of zeros\n');

else 

    ZeroNum=0;

end    



%生成信号，signal是原信号。signal为采样信号。

fs=1/T;

t=0:0.00001:T*(N+ZeroNum-1);

signal=sin(2*pi*f*t); 

t2=0:T:T*(N+ZeroNum-1);

signal2=sin(2*pi*f*t2);

if (flag)

    signal2=[signal2 zeros(1, ZeroNum)];

end



%画出原信号及采样信号。

figure;

subplot(2,1,1);

plot(t,signal);

xlabel('Time(s)');

ylabel('Amplitude(volt)');

title('Singnal');

hold on;

subplot(2,1,1);

stem(t2,signal2,'r');

axis([0 T*(N+ZeroNum) -1 1]);



%作FFT变换，计算其幅值，归一化处理，并画出频谱。

Y = fft(signal2,N);

Pyy = Y.* conj(Y) ;

Pyy=(Pyy/sum(Pyy))*2;

f=0:fs/(N-1):fs/2;4

subplot(2,1,2);

bar(f,Pyy(1:N/2));

xlabel('Frequency(Hz)');

ylabel('Amplitude');

title('Frequency compnents of signal');

axis([0 fs/2 0 ceil(max(Pyy))])

grid on;

1万 · 2015-12-29 22:48

复制

FFT与反FFT的源码：



/*

Parameter:

f:        point a source sequance f(x)

F:        point a target sequance F(x) by FFT procedure

Function:

the procedure through FFT transform f(x) to F(x)

*/

void CJPEG::FFT(const complex<double> *f, complex<double> *F, int N)

{

    

    //将时域复制至频域

    memcpy(F,f,sizeof(complex<double>) * N);



    //需要多少级蝶形运算

    int r = log(N +1) / log(2);



    //进行倒位序运算

    BitReOrder(F,r);



    //计算Wr因子, 加权系数 

    complex<double> *W = new complex<double>[N/2];

    double angle;    

    for(int i = 0; i < N / 2; i++) 

    { 

        angle = - i * PI * 2 / N; 

        W[i] = complex<double> (cos(angle), sin(angle)); 

        //W[i] = complex<double>(0, exp(-i*2*PI/N));

    } 



    //采用蝶形算法进行快速傅立叶变换

    int DFTn;//第DFTn级

    int k;//分级有多少个蝶形运算

    int d;//蝶形运算的偏移

    int p;//index

    

    //complex<double> *X,*X1,*X2;

    //X1 = new complex<double>[N]; 

    //X2 = new complex<double>[N]; 

    complex<double> X1,X2;

    for( DFTn = 0; DFTn < r; DFTn++){//第几级蝶形运算

        for( k= 0; k < (1 << (r - DFTn - 1)); ++k){//当前列所有的蝶形进行运算

            p = 2 * k *  (1 << DFTn);

            for ( d = 0; d < (1 << DFTn); ++d){//相邻蝶形运算

                //p = k * (1 << (k + 1)) + d;

                X1 = F[p+d];

                X2 = F[p+d+(1 << DFTn)];

                F[p+d] = X1 + X2 * W[d*(1<<(r - DFTn - 1))];

                F[p+d+(1 << DFTn)] = X1 - X2 * W[d*(1<<(r - DFTn - 1))];

            }

        }

    }

    //BitReOrder(F,r);

    for(i = 0; i < N; ++i)

        F[i] /= N;

    delete[] W;



    return;

}







void CJPEG::IFFT(const complex<double> *F, complex<double> *f, int N)

{

    int i;

    complex<double>*T = new complex<double>[N];

    memcpy(T,F,sizeof(complex<double>) * N);

    for(i = 0; i < N; ++i){

        //取共轭

        T[i] = complex<double>(T[i].real(), -T[i].imag());        

    }

    CJPEG::FFT(T,f,N);

    for(i = 0; i < N; ++i){//取共轭，并乘以N

        f[i] = complex<double>(f[i].real(), -f[i].imag());    

        f[i] *= complex<double>(N,0) ;

    }

    delete[] T;

}



int CJPEG::ReadFraqData(const char *lpszFileName, int N, double *pData)

{

    return 0;

}



int CJPEG::ReadTimeData(const char *lpszFileName, int N, double *pData)

{

    if(pData == NULL) return -1;

    CTextFile ctf(lpszFileName);

    char sLine[64];

    double dData;

    int i = 0;

    while(!ctf.Eof() && i < N){

        ctf.GetLine(sLine, 64);

        dData = atof(sLine);

        pData[i++] = dData;

    }

    return i;

}



/////////////////////////////////////////// 

// Function name : BitReverse 

// Description : 二进制倒序操作 

// Return type : int 

// Argument : int src 待倒读的数 

// Argument : int size 二进制位数 

int CJPEG::BitReverse(int src, int size) 

{ 

    int tmp = src; 

    int des = 0; 

    for (int i=size-1; i>=0; i--) 

    { 

        des = ((tmp & 0x1) << i) | des; 

        tmp = tmp >> 1; 

    } 

    return des; 

} 



/*

Parameter:

sequ:    point a sequance that f(x)

N:        f(x), the x maxinum

Function:

the procedure change the sequance by this:

    // 101 -> 101

    // 110 -> 011

    // 1010 -> 0101

f(5) -> f(5)

f(6) -> f(3)

f(10) -> f(5)

and so on.

*/

void CJPEG::BitReOrder(complex<double> *sequ, int r)

{

    complex<double> dTemp;

    int x;

    for(int i = 0; i <(1 << r); ++i){

        x = BitReverse(i,r);

        if (x > i){

            dTemp = sequ[i];

            sequ[i] = sequ[x];

            sequ[x] = dTemp;

        }

    }    

}

1万 · 2016-1-8 20:41

#include "math.h"
#define FFT_N 256             //----N点FFT变换
#define FFT_M 8             //----N等于2的M次方

u16 weizhi[FFT_N];
u16 wei1[FFT_N];
float temdat1[FFT_N];
float temdat2[FFT_N];
float tem_shi[FFT_N];
float tem_xu[FFT_N];
float out_shi[FFT_N];
float out_xu[FFT_N];
float Pix=3.1415926;
float P=0;
float v1,v2;
float v3,v4,v5,v6;

//----------------------------------------------------------------------------------
//--- 函数
//f=(double)i/100; //-- i的取值为 0-255，步进为1
//ddt=sin(2*pi*20*f);
//--正弦函数，频率20Hz//--以fs=100Hz的采样频率采样，采样N=256个点，进行256点FFT运算
//----在数据点位x=51时--频率量的值达到最大点此时的频率量为 f=(x*fs)/(N-1)=20Hz
//----------------------------------------------------------------------------------
float ddt[FFT_N]=
{

0.000000,0.951057,0.587785,-0.587785,-0.951057,-0.000000,0.951056,0.587785,
-0.587785,-0.951057,-0.000000,0.951056,0.587785,-0.587785,-0.951057,-0.000000,
0.951056,0.587786,-0.587785,-0.951057,-0.000000,0.951056,0.587786,-0.587785,
-0.951057,-0.000001,0.951056,0.587786,-0.587785,-0.951057,-0.000001,0.951056,
0.587786,-0.587785,-0.951057,-0.000001,0.951056,0.587786,-0.587785,-0.951057,
-0.000001,0.951056,0.587786,-0.587785,-0.951057,-0.000001,0.951056,0.587786,
-0.587784,-0.951057,-0.000001,0.951056,0.587786,-0.587784,-0.951057,-0.000001,
0.951056,0.587786,-0.587784,-0.951057,-0.000001,0.951056,0.587786,-0.587784,
-0.951057,-0.000001,0.951056,0.587786,-0.587784,-0.951057,-0.000002,0.951056,
0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587784,-0.951057,
-0.000002,0.951056,0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,
-0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587784,-0.951057,-0.000002,
0.951056,0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587783,
-0.951057,-0.000002,0.951056,0.587787,-0.587783,-0.951057,-0.000002,0.951056,
0.587787,-0.587783,-0.951057,-0.000002,0.951056,0.587787,-0.587783,-0.951057,
-0.000003,0.951056,0.587787,-0.587783,-0.951057,-0.000003,0.951056,0.587787,
-0.587783,-0.951057,-0.000003,0.951056,0.587788,-0.587783,-0.951057,-0.000003,
0.951056,0.587788,-0.587783,-0.951057,-0.000003,0.951056,0.587788,-0.587783,
-0.951057,-0.000003,0.951056,0.587788,-0.587783,-0.951058,-0.000003,0.951056,
0.587788,-0.587783,-0.951058,-0.000003,0.951055,0.587788,-0.587783,-0.951058,
-0.000003,0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,0.587788,
-0.587782,-0.951058,-0.000004,0.951055,0.587788,-0.587782,-0.951058,-0.000004,
0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,0.587788,-0.587782,
-0.951058,-0.000004,0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,
0.587789,-0.587782,-0.951058,-0.000004,0.951055,0.587789,-0.587782,-0.951058,
-0.000004,0.951055,0.587789,-0.587782,-0.951058,-0.000004,0.951055,0.587789,
-0.587782,-0.951058,-0.000005,0.951055,0.587789,-0.587782,-0.951058,-0.000005,
0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587789,-0.587781,
-0.951058,-0.000005,0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,
0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587789,-0.587781,-0.951058,
-0.000005,0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587790,
-0.587781,-0.951058,-0.000005,0.951055,0.587790,-0.587781,-0.951058,-0.000005,

};

//------------------------------------------------------------------------------
//-- 复数求模  ----------------------------------------------------------------
//------------------------------------------------------------------------------
float DAT_sqrt(float A1,float A2)
{
u8 i=0;
float x;
      float y;
      x=A1*A1;
y=A2*A2;
      y=x+y;
if(y<1) x=1;
else  x=y/2;
for(i=0;i<15;i++)
{
   if(x*x<=y)
                              break;
   else
      x=((y/x)+x)/2; //牛顿迭代法
}
return x;
}

//-----------------------------------------------------------------------------
//---快速傅里叶变换 N个点变换--M层蝶形计算
//-----------------------------------------------------------------------------
void FFT(float in_dat[],float *jieguo)
{
u16 L,B;
u16 i=0,j=0;
u16 j1,j2=0;
      u16 N=FFT_N;
      u16 M=FFT_M;

//-----------------------------位置置换 b2b1b0---b0b1b2--做M次左移
for(i=0;i<N;i++)
{
   weizhi[i]=i;
   wei1[i]=0;
}
for(i=0;i<N;i++)
{
   for(j=0;j<M;j++)
   {
      wei1[i]=wei1[i]<<1;
      if(weizhi[i]&(0x01<<j))
      wei1[i]|=0x01;
      else
      wei1[i]&=0xfffe;
   }
   weizhi[i]=wei1[i];
}
for(i=0;i<N;i++)
temdat1[i]=in_dat[i];
for(i=0;i<N;i++)
temdat2[i]=temdat1[weizhi[i]];
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//----变换位置的数据存在 temdat2[]中
//---- 欧拉公式 e^(-jx)=cos(x)-jsin(x);
//------蝶形计算----------------------------------------------------------------
//---B表示蝶形中两点的距离---L表示第L级蝶形运算--j表示的j行的计算
//------------------------------------------------------------------------------
  for(i=0;i<N;i++)
  {
out_shi[i]=temdat2[i];
out_xu[i]=0.0;
  }
  B=1;

  for(L=1;L<=M;L++)
  {
for(j=0;j<N;j++)
{
   j1=1<<L;
   j2=1<<(L-1);
   if((j%j1)<j2)
   {
   P=j*N;
   for(i=0;i<(L-1);i++)
   P=P/2;
   P=P*Pix;
   P=P/N;
   v3=cos(P);
   v4=0-sin(P);
   v5=(out_shi[j+B]*v3);
   v6=(out_xu[j+B]*v4);
   v1=v5-v6;
   v5=(out_shi[j+B]*v4);
   v6=(out_xu[j+B]*v3);
   v2=v5+v6;
   tem_shi[j]  =out_shi[j]+v1;       //--计算实部
   tem_xu[j] =out_xu[j]+v2;       //--计算虚部
   tem_shi[j+B]=out_shi[j]-v1;       //--计算实部
   tem_xu[j+B] =out_xu[j]-v2;       //--计算虚部
   out_shi[j]=tem_shi[j];
   out_xu[j] =tem_xu[j];
   out_shi[j+B]=tem_shi[j+B];
   out_xu[j+B] =tem_xu[j+B];
   }
}
B=B*2;
  }
//-------------------------------------------------------
//---FFT 变换后的复数取模--------------------------------
//-------------------------------------------------------
      for(i=0;i<N;i++)
      {
            *jieguo=DAT_sqrt(out_shi[i],out_xu[i]);
            jieguo++;
      }
//-------------------------------------------------------
}

//----------------------------------------------
//--变换结果如下
//----------------------------------------------
void shujuchuli(void)
{

      FFT(ddt,ddt);

}

MATLAB软件计算程序

%---------------------------------------------------------
clf;
fs=100;N=256; %采样频率和数据点数

n=0:N-1;t=n/fs; %时间序列

x=sin(2*pi*20*t); %信号

y=fft(x,N); %对信号进行快速Fourier变换
mag=abs(y);    %求得Fourier变换后的振幅

f=n*fs/N; %频率序列

subplot(2,2,3),plot(n,x); %绘出随频率变化的振幅
xlabel('t/Hz');
ylabel('振幅');title('N=256');grid on;

subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=256');grid on;

subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=256');grid on;

%---------------------------------------------------------