- #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;
- %---------------------------------------------------------