[应用相关] STM32--任意点数FFT变换程序

[复制链接]
 楼主| 734774645 发表于 2015-12-29 22:27 | 显示全部楼层 |阅读模式
  1. #include "math.h"
  2. #define FFT_N 256                //----N点FFT变换
  3. #define FFT_M 8                //----N等于2的M次方

  4. u16 weizhi[FFT_N];
  5. u16 wei1[FFT_N];
  6. float temdat1[FFT_N];
  7. float temdat2[FFT_N];
  8. float tem_shi[FFT_N];
  9. float tem_xu[FFT_N];  
  10. float out_shi[FFT_N];
  11. float out_xu[FFT_N];  
  12. float Pix=3.1415926;
  13. float P=0;
  14. float v1,v2;
  15. float v3,v4,v5,v6;



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

  25. 0.000000,0.951057,0.587785,-0.587785,-0.951057,-0.000000,0.951056,0.587785,
  26. -0.587785,-0.951057,-0.000000,0.951056,0.587785,-0.587785,-0.951057,-0.000000,
  27. 0.951056,0.587786,-0.587785,-0.951057,-0.000000,0.951056,0.587786,-0.587785,
  28. -0.951057,-0.000001,0.951056,0.587786,-0.587785,-0.951057,-0.000001,0.951056,
  29. 0.587786,-0.587785,-0.951057,-0.000001,0.951056,0.587786,-0.587785,-0.951057,
  30. -0.000001,0.951056,0.587786,-0.587785,-0.951057,-0.000001,0.951056,0.587786,
  31. -0.587784,-0.951057,-0.000001,0.951056,0.587786,-0.587784,-0.951057,-0.000001,
  32. 0.951056,0.587786,-0.587784,-0.951057,-0.000001,0.951056,0.587786,-0.587784,
  33. -0.951057,-0.000001,0.951056,0.587786,-0.587784,-0.951057,-0.000002,0.951056,
  34. 0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587784,-0.951057,
  35. -0.000002,0.951056,0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,
  36. -0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587784,-0.951057,-0.000002,
  37. 0.951056,0.587787,-0.587784,-0.951057,-0.000002,0.951056,0.587787,-0.587783,
  38. -0.951057,-0.000002,0.951056,0.587787,-0.587783,-0.951057,-0.000002,0.951056,
  39. 0.587787,-0.587783,-0.951057,-0.000002,0.951056,0.587787,-0.587783,-0.951057,
  40. -0.000003,0.951056,0.587787,-0.587783,-0.951057,-0.000003,0.951056,0.587787,
  41. -0.587783,-0.951057,-0.000003,0.951056,0.587788,-0.587783,-0.951057,-0.000003,
  42. 0.951056,0.587788,-0.587783,-0.951057,-0.000003,0.951056,0.587788,-0.587783,
  43. -0.951057,-0.000003,0.951056,0.587788,-0.587783,-0.951058,-0.000003,0.951056,
  44. 0.587788,-0.587783,-0.951058,-0.000003,0.951055,0.587788,-0.587783,-0.951058,
  45. -0.000003,0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,0.587788,
  46. -0.587782,-0.951058,-0.000004,0.951055,0.587788,-0.587782,-0.951058,-0.000004,
  47. 0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,0.587788,-0.587782,
  48. -0.951058,-0.000004,0.951055,0.587788,-0.587782,-0.951058,-0.000004,0.951055,
  49. 0.587789,-0.587782,-0.951058,-0.000004,0.951055,0.587789,-0.587782,-0.951058,
  50. -0.000004,0.951055,0.587789,-0.587782,-0.951058,-0.000004,0.951055,0.587789,
  51. -0.587782,-0.951058,-0.000005,0.951055,0.587789,-0.587782,-0.951058,-0.000005,
  52. 0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587789,-0.587781,
  53. -0.951058,-0.000005,0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,
  54. 0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587789,-0.587781,-0.951058,
  55. -0.000005,0.951055,0.587789,-0.587781,-0.951058,-0.000005,0.951055,0.587790,
  56. -0.587781,-0.951058,-0.000005,0.951055,0.587790,-0.587781,-0.951058,-0.000005,


  57. };


  58. //------------------------------------------------------------------------------
  59. //-- 复数 求模  ----------------------------------------------------------------
  60. //------------------------------------------------------------------------------
  61. float DAT_sqrt(float A1,float A2)
  62. {
  63.    u8 i=0;
  64.    float x;
  65.          float y;
  66.          x=A1*A1;
  67.    y=A2*A2;
  68.          y=x+y;
  69.    if(y<1) x=1;
  70.    else  x=y/2;  
  71.    for(i=0;i<15;i++)
  72.    {   
  73.      if(x*x<=y)
  74.                                 break;
  75.      else
  76.         x=((y/x)+x)/2; //牛顿迭代法  
  77.    }
  78.    return x;
  79. }


  80. //-----------------------------------------------------------------------------
  81. //---快速傅里叶变换 N个点变换--M层 蝶形计算
  82. //-----------------------------------------------------------------------------
  83. void FFT(float in_dat[],float *jieguo)
  84. {
  85.    u16 L,B;
  86.    u16 i=0,j=0;
  87.    u16 j1,j2=0;
  88.          u16 N=FFT_N;
  89.          u16 M=FFT_M;
  90.         
  91. //-----------------------------位置置换 b2b1b0---b0b1b2--做M次左移
  92.    for(i=0;i<N;i++)
  93.    {
  94.      weizhi[i]=i;
  95.      wei1[i]=0;
  96.    }
  97.    for(i=0;i<N;i++)
  98.    {
  99.       for(j=0;j<M;j++)
  100.       {
  101.         wei1[i]=wei1[i]<<1;
  102.         if(weizhi[i]&(0x01<<j))
  103.         wei1[i]|=0x01;
  104.         else
  105.         wei1[i]&=0xfffe;              
  106.       }
  107.       weizhi[i]=wei1[i];                     
  108.    }
  109.    for(i=0;i<N;i++)
  110.    temdat1[i]=in_dat[i];
  111.    for(i=0;i<N;i++)
  112.    temdat2[i]=temdat1[weizhi[i]];
  113. //------------------------------------------------------------------------------   
  114. //------------------------------------------------------------------------------
  115. //----变换位置的数据存在 temdat2[]中
  116. //---- 欧拉公式 e^(-jx)=cos(x)-jsin(x);     
  117. //------蝶形计算----------------------------------------------------------------
  118. //---B表示蝶形中两点的距离---L表示第L级蝶形运算--j表示的j行的计算
  119. //------------------------------------------------------------------------------
  120.   for(i=0;i<N;i++)
  121.   {
  122.    out_shi[i]=temdat2[i];
  123.    out_xu[i]=0.0;
  124.   }
  125.   B=1;
  126.                                 
  127.   for(L=1;L<=M;L++)
  128.   {
  129.    for(j=0;j<N;j++)
  130.    {
  131.      j1=1<<L;
  132.      j2=1<<(L-1);
  133.      if((j%j1)<j2)
  134.      {
  135.       P=j*N;
  136.       for(i=0;i<(L-1);i++)
  137.       P=P/2;
  138.       P=P*Pix;
  139.       P=P/N;
  140.       v3=cos(P);
  141.       v4=0-sin(P);
  142.       v5=(out_shi[j+B]*v3);
  143.       v6=(out_xu[j+B]*v4);   
  144.       v1=v5-v6;        
  145.       v5=(out_shi[j+B]*v4);
  146.       v6=(out_xu[j+B]*v3);
  147.       v2=v5+v6;                              
  148.       tem_shi[j]  =out_shi[j]+v1;         //--计算实部                  
  149.       tem_xu[j]   =out_xu[j]+v2;          //--计算虚部  
  150.       tem_shi[j+B]=out_shi[j]-v1;         //--计算实部      
  151.       tem_xu[j+B] =out_xu[j]-v2;          //--计算虚部  
  152.       out_shi[j]=tem_shi[j];
  153.       out_xu[j] =tem_xu[j];
  154.       out_shi[j+B]=tem_shi[j+B];
  155.       out_xu[j+B] =tem_xu[j+B];                                    
  156.      }            
  157.    }
  158.    B=B*2;                 
  159.   }
  160. //-------------------------------------------------------
  161. //---FFT 变换后的复数取模--------------------------------
  162. //-------------------------------------------------------
  163.         for(i=0;i<N;i++)
  164.         {
  165.                 *jieguo=DAT_sqrt(out_shi[i],out_xu[i]);
  166.                 jieguo++;
  167.         }
  168. //-------------------------------------------------------   
  169. }



  170. //----------------------------------------------
  171. //--变换结果如下
  172. //----------------------------------------------
  173. void shujuchuli(void)
  174. {

  175.         FFT(ddt,ddt);


  176. }




  177. MATLAB软件计算程序

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

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

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

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

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

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


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


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

  195. %---------------------------------------------------------




598330983 发表于 2015-12-29 22:47 | 显示全部楼层
一,如果对信号进行同样点数N的FFT变换,采样频率fs越高,则可以分析越高频的信号;与此同时,采样频率越低,对于低频信号的频谱分辨率则越好。

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

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

四,如果采样频率fs小于2倍信号频率2*fs(奈圭斯特定理),则频谱分析结果会出错。
598330983 发表于 2015-12-29 22:47 | 显示全部楼层
  1. %清除命令窗口及变量
  2. clc;
  3. clear all;

  4. %输入f、N、T、是否补零(补几个零)
  5. f=input('Input frequency of the signal: f\n');
  6. N=input('Input number of pointsl: N\n');
  7. T=input('Input sampling time: T\n');
  8. flag=input('Add zero too sampling signal or not? yes=1 no=0\n');
  9. if(flag)
  10.     ZeroNum=input('Input nmber of zeros\n');
  11. else
  12.     ZeroNum=0;
  13. end   

  14. %生成信号,signal是原信号。signal为采样信号。
  15. fs=1/T;
  16. t=0:0.00001:T*(N+ZeroNum-1);
  17. signal=sin(2*pi*f*t);
  18. t2=0:T:T*(N+ZeroNum-1);
  19. signal2=sin(2*pi*f*t2);
  20. if (flag)
  21.     signal2=[signal2 zeros(1, ZeroNum)];
  22. end

  23. %画出原信号及采样信号。
  24. figure;
  25. subplot(2,1,1);
  26. plot(t,signal);
  27. xlabel('Time(s)');
  28. ylabel('Amplitude(volt)');
  29. title('Singnal');
  30. hold on;
  31. subplot(2,1,1);
  32. stem(t2,signal2,'r');
  33. axis([0 T*(N+ZeroNum) -1 1]);

  34. %作FFT变换,计算其幅值,归一化处理,并画出频谱。
  35. Y = fft(signal2,N);
  36. Pyy = Y.* conj(Y) ;
  37. Pyy=(Pyy/sum(Pyy))*2;
  38. f=0:fs/(N-1):fs/2;4
  39. subplot(2,1,2);
  40. bar(f,Pyy(1:N/2));
  41. xlabel('Frequency(Hz)');
  42. ylabel('Amplitude');
  43. title('Frequency compnents of signal');
  44. axis([0 fs/2 0 ceil(max(Pyy))])
  45. grid on;


598330983 发表于 2015-12-29 22:48 | 显示全部楼层
  1. FFT与反FFT的源码:

  2. /*
  3. Parameter:
  4. f:        point a source sequance f(x)
  5. F:        point a target sequance F(x) by FFT procedure
  6. Function:
  7. the procedure through FFT transform f(x) to F(x)
  8. */
  9. void CJPEG::FFT(const complex<double> *f, complex<double> *F, int N)
  10. {
  11.    
  12.     //将时域复制至频域
  13.     memcpy(F,f,sizeof(complex<double>) * N);

  14.     //需要多少级蝶形运算
  15.     int r = log(N +1) / log(2);

  16.     //进行倒位序运算
  17.     BitReOrder(F,r);

  18.     //计算Wr因子, 加权系数
  19.     complex<double> *W = new complex<double>[N/2];
  20.     double angle;   
  21.     for(int i = 0; i < N / 2; i++)
  22.     {
  23.         angle = - i * PI * 2 / N;
  24.         W[i] = complex<double> (cos(angle), sin(angle));
  25.         //W[i] = complex<double>(0, exp(-i*2*PI/N));
  26.     }

  27.     //采用蝶形算法进行快速傅立叶变换
  28.     int DFTn;//第DFTn级
  29.     int k;//分级有多少个蝶形运算
  30.     int d;//蝶形运算的偏移
  31.     int p;//index
  32.    
  33.     //complex<double> *X,*X1,*X2;
  34.     //X1 = new complex<double>[N];
  35.     //X2 = new complex<double>[N];
  36.     complex<double> X1,X2;
  37.     for( DFTn = 0; DFTn < r; DFTn++){//第几级蝶形运算
  38.         for( k= 0; k < (1 << (r - DFTn - 1)); ++k){//当前列所有的蝶形进行运算
  39.             p = 2 * k *  (1 << DFTn);
  40.             for ( d = 0; d < (1 << DFTn); ++d){//相邻蝶形运算
  41.                 //p = k * (1 << (k + 1)) + d;
  42.                 X1 = F[p+d];
  43.                 X2 = F[p+d+(1 << DFTn)];
  44.                 F[p+d] = X1 + X2 * W[d*(1<<(r - DFTn - 1))];
  45.                 F[p+d+(1 << DFTn)] = X1 - X2 * W[d*(1<<(r - DFTn - 1))];
  46.             }
  47.         }
  48.     }
  49.     //BitReOrder(F,r);
  50.     for(i = 0; i < N; ++i)
  51.         F[i] /= N;
  52.     delete[] W;

  53.     return;
  54. }



  55. void CJPEG::IFFT(const complex<double> *F, complex<double> *f, int N)
  56. {
  57.     int i;
  58.     complex<double>*T = new complex<double>[N];
  59.     memcpy(T,F,sizeof(complex<double>) * N);
  60.     for(i = 0; i < N; ++i){
  61.         //取共轭
  62.         T[i] = complex<double>(T[i].real(), -T[i].imag());        
  63.     }
  64.     CJPEG::FFT(T,f,N);
  65.     for(i = 0; i < N; ++i){//取共轭,并乘以N
  66.         f[i] = complex<double>(f[i].real(), -f[i].imag());   
  67.         f[i] *= complex<double>(N,0) ;
  68.     }
  69.     delete[] T;
  70. }

  71. int CJPEG::ReadFraqData(const char *lpszFileName, int N, double *pData)
  72. {
  73.     return 0;
  74. }

  75. int CJPEG::ReadTimeData(const char *lpszFileName, int N, double *pData)
  76. {
  77.     if(pData == NULL) return -1;
  78.     CTextFile ctf(lpszFileName);
  79.     char sLine[64];
  80.     double dData;
  81.     int i = 0;
  82.     while(!ctf.Eof() && i < N){
  83.         ctf.GetLine(sLine, 64);
  84.         dData = atof(sLine);
  85.         pData[i++] = dData;
  86.     }
  87.     return i;
  88. }

  89. ///////////////////////////////////////////
  90. // Function name : BitReverse
  91. // Description : 二进制倒序操作
  92. // Return type : int
  93. // Argument : int src 待倒读的数
  94. // Argument : int size 二进制位数
  95. int CJPEG::BitReverse(int src, int size)
  96. {
  97.     int tmp = src;
  98.     int des = 0;
  99.     for (int i=size-1; i>=0; i--)
  100.     {
  101.         des = ((tmp & 0x1) << i) | des;
  102.         tmp = tmp >> 1;
  103.     }
  104.     return des;
  105. }

  106. /*
  107. Parameter:
  108. sequ:    point a sequance that f(x)
  109. N:        f(x), the x maxinum
  110. Function:
  111. the procedure change the sequance by this:
  112.     // 101 -> 101
  113.     // 110 -> 011
  114.     // 1010 -> 0101
  115. f(5) -> f(5)
  116. f(6) -> f(3)
  117. f(10) -> f(5)
  118. and so on.
  119. */
  120. void CJPEG::BitReOrder(complex<double> *sequ, int r)
  121. {
  122.     complex<double> dTemp;
  123.     int x;
  124.     for(int i = 0; i <(1 << r); ++i){
  125.         x = BitReverse(i,r);
  126.         if (x > i){
  127.             dTemp = sequ[i];
  128.             sequ[i] = sequ[x];
  129.             sequ[x] = dTemp;
  130.         }
  131.     }   
  132. }


 楼主| 734774645 发表于 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;

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


mintspring 发表于 2016-1-8 22:35 | 显示全部楼层
对于同样的采样率fs,提高FFT的点数N,则可提高频谱的分辨率。
网络孤客 发表于 2016-1-10 14:29 | 显示全部楼层
谢谢楼主,慢慢研究。
gejigeji521 发表于 2016-1-31 11:04 | 显示全部楼层
对于同样的采样率fs,提高FFT的点数N,则可提高频谱的分辨率。
ham56 发表于 2016-11-18 18:47 来自手机 | 显示全部楼层
谢谢楼主,慢慢理解,学习
sxs20122730 发表于 2018-5-30 16:48 | 显示全部楼层
我把程序在keil上仿真,256点没问题,但是512点的FFT就一直跑不通,请问之前你有试过512点的FFT吗?期待你的回复
akywong 发表于 2021-3-19 18:51 | 显示全部楼层
sxs20122730 发表于 2018-5-30 16:48
我把程序在keil上仿真,256点没问题,但是512点的FFT就一直跑不通,请问之前你有试过512点的FFT吗?期待你 ...

可以,跑8K个点都可以的
您需要登录后才可以回帖 登录 | 注册

本版积分规则

210

主题

3585

帖子

15

粉丝
快速回复 在线客服 返回列表 返回顶部