一、FFT是什么，能干啥？
% n' ^  o( ^" x6 V$ K2 p7 [FFT(fast Fourier transform)快速傅里叶变换，把时域信号变换到频域，至于换到频域后能做啥，那就多了,你能搜到这篇文章，说明你已经知道他是干啥的了，，，，
! _! z- B+ B; G; v
  k( Q8 Z! c! ^
8 O; K! W" _  U0 S  c6 B2 C# D二、先用matlab来认识一下FFT
) i5 {2 S) R8 D$ w: k. x1.上代码
代码如下：

1. fs = 1024;           %采样频率
   8 \9 ]9 B$ J9 l8 G. H
2. NPT = 1024;          %样点个数
   ( Q3 j5 K3 n$ }% E0 K
3. t=0:1/fs:1;          %采样步长,模拟单片机采样频率
   2 B* C. H) D% J& O0 P7 E
4. y= 2+200*cos(2*pi*15*t)+40*cos(2*pi*100*t+pi*90/180)+15*cos(2*pi*250*t+pi*90/180);
   % x. v0 ]0 O8 ^
5. plot(t,y);                     %画出未变换前的时域信号波形
   
6. fbl=fs/NPT;          %分辨率
   7 {7 V  S. D4 n: c; ^
7. f=(0:NPT)*fbl;       %每个点所代表的频率
   
8. Y=fft(y);                     %进行FFT变换，matlab中fft函数，输入多少采样数据，就输出多少幅频值(是个复数)
   , F0 U: k! c: u$ y: ]+ k# Z* J! E$ i
9.                                      %所以length(y)=length(Y);
   
10. figure(2)
    
11. plot(f(1:NPT/2),abs(Y(1:NPT/2)));
    + g& x6 L3 P* g; Q* w
12.                                      %以每个点所代表的频率画出其对应的赋值
    
13. 
    3 G, e# _1 p& n% d5 W8 j7 D

复制代码

2.代码关键字讲解
思想都是一样的，matab知道了，在32里面也差不多
, l' O9 ~/ \# d7 }  z* Q3 C3 i
采样频率： 用单片机ADC采样，结果是离散的数字信号，参考采样定理，采样频率至少大于信号频率的两倍，在这里我们采样频率1024,远大于基波频率2，选取1024也是为了方便自己，可以更直观的观察输出（2的整数倍都可）。
分辨率： 输出时，最小频率间隔(Fs/N)，第n个点所代表的频率即为Fn=(n-1)*Fs/N,在这个程序中，采样频率是1024，采样点数是1024，所以，分辨率为1Hz
5 I6 @+ o: l! C8 F5 Q5 d输出赋值 fft输出就是复数，取模值即可
" I! {8 f! H+ ^8 Z/ t! u0 B8 y& |
  h- _2 V; w; @$ s1 X3 \1 Q
: q' B9 }+ y4 f8 N; P3.程序结果
+ q" w# g. \- m2 o未变换前的时域信号波形图：

. b, _6 R0 z$ W

4 D# r# g! p6 h3 @转换后输出的幅频函数波形图：
% x: g6 V* o  a! r- _( t9 }
, W7 L+ g7 F: e7 o

" i0 S( F0 z% l& t: p$ A由图可看出输入信号大概有三个频率，15，100，250.与输入信号相符
二、淦! 上32
; B# ~$ i# g2 C7 A) m' K& eSTM32官方给了三个函数库文件，分别是：
+ ]8 @% F2 @/ L& _cr4_fft_64_stm32.s
cr4_fft_256_stm32.s
, E! Z) R& O' ^: h+ I/ jcr4_fft_1024_stm32.s
5 t. C1 C8 Z- q1 s* d; k/ Z（这些都是汇编语言写的，好在我们不需要看懂，只需要会调用就好。）
- m/ \) x4 a' W3 O分别对应这采样点数为64，256，1024。库函数的获取，可以在官网上获取我们的目标是计算周期为1ms波形的THD值！，使用采样点数为64位的函数

" }& I" W& o0 h. }0 g1、代码移植
9 O$ F2 e; h' P( P9 R新建文件夹，把
1 a3 A, u1 b& f2 P) o- O& Ccr4_fft_64_stm32.s
, n+ A" j' A) acr4_fft_256_stm32.s
cr4_fft_1024_stm32.s
和stm32_dsp.h
文件添加到工程里（就像添加.C文件一样），注：stm32_dsp.h文件里第27行需要根据自己单片机型号来修改
8 a1 G. B3 R; k- W- e% u
# m2 C2 H& Y4 l, b) Z; ~; D
# u! |, z: _! L) i

2，ADC采集
7 e, J7 V$ a- n1 ?/ P: JADC采集没有什么好说的，我想让大家注意一点就是采样周期，在这里，为了方便使用，我们设ADC采样频率为64KHz，注意，这里直接调用正点原子的Get_Adc();函数的话，他的执行周期很长，需要改短（自己跳进Get_Adc(）函数里，看函数定义，再改），要不然，在定时器中，ADC采集时间会超过中断事件，导致采样周期不是64KHz,最好在中断中设置一个LED灯反转的程序，检测ADC采样周期，确保ADC采集频率稳定

+ n! ]! V$ Q! h) f. ]
+ R# J: O1 f1 O4 {" W: h3 u, Q  K8 E3、上函数
: F% G1 x2 l3 A在上函数之前，我们把函数的每一个输入输出变量的含义说一下
: }$ m, L$ ]% m! @8 l在我的程序里是这样调用的：cr4_fft_64_stm32(output,input,NPT);
她的输入变量是input：是一个int32_t类型的数组，高十六位是实部，第十六位是虚部，显然，我们ADC采集只能采集实数。所以我对采集到的数据进行这样的处理：
4 G) ^4 G4 E+ E, y% z; X

1.   ADC_Value[i]=Get_Adc(ADC_Channel_5);
   + w8 N! \& I& p
2.          input[i] = ((s16)ADC_Value[i])<<16;

复制代码

至于虚部，就不用管啦，移位后默认是0
* K+ E. U; G( u; i3 {! p  \: ~他的另一个输入是NPT ，采样点数。最好是定义成宏定义，因为要用到他定义的数组太多，这里是64
他的输出是output,也是个复数，高十六位代表实部，低十六位代表虚部，他们的模值就是我们要的幅值
呐，官方给了专门求赋值的函数
% p2 x2 _! q) v# M
! M( L' o2 ?7 o% _8 a

1. void dsp_asm_powerMag(int32_t *IBUFOUT)
   ; G) c5 H0 }- {) j
2. {
   9 o- k9 P5 z2 r
3.         s16 lx,ly;
   , ~1 ]. l: q; J. w3 Y2 D* H
4.         u32 i;
   5 e4 g1 j* l" ~0 g+ Y5 A2 U
5.         for(i=0;i<NPT/2;i++)//由于FFT的频域结果是关于奈奎斯特频率对称的，所以只计算一半
   
6.         {
   ; m6 I( i: u9 S
7.                 lx = (IBUFOUT[i] << 16)>>16;      //取低十六位，虚部
   : ^3 O, ^9 ?7 _  [
8.                 ly = (IBUFOUT[i] >> 16);          //取高十六位，实部
   2 m1 M% A2 ^4 A2 R2 j$ {4 u
9.                 float X = NPT*((float)lx)/32768;
   
10.                 float Y = NPT*((float)ly)/32768;
    4 ?0 r% ]+ r# s4 m* ^6 {' h0 H
11.                 float Mag = sqrt(X*X+Y*Y)/NPT;
    
12.                 OUTPUT_MAG[i] = (u32)(Mag*65536);
      |& o4 I( `& N0 W
13.         }//这些就是计算振幅IBUFOUT[i] = sqrt(lx*lx+ly*ly)*2/NPT
    3 p1 @6 T& V% I6 {% b; J
14.          //之所以先除以32768后乘以65536是为了符合浮点数的计算规律
    $ e& R+ W6 K* A; l/ R& H' {# C
15.         OUTPUT_MAG[0] = OUTPUT_MAG[0]/2;//这个是直流分量，不需要乘以2
    
16. }
    3 L2 [" {: D8 A% ~0 z0 r# I

复制代码

直接用就很OK了

1. dsp_asm_powerMag(output);               
   

复制代码

& J- n) r3 ]/ ?) t& `9 d那么我们调用了函数，也转换成幅值了，所以，我该怎么用这些输出的数呢。他输出的是什么***，接下来就需要对应起来，
我们需要知道的是 每个频率点对应的幅值，现在有幅值，怎么把他对应的点找到，这样我们就大功告成了
6 V; Z  i6 _/ D0 f1 d% a4 N" ?6 c还记得前面讲的分辨率嘛，对，我们采样点是64个，采样频率也是64KHz,那么分辨率就是fs/NPT即为1KHz，也就是说输出的output每隔一个对应着一个相差位1Hz的频率，那么output[1]就对应着1KHz，output[2]就对应着2KHz，output[3]就对应着3KHz，，，，我们只取前NPT/2个，不要问为什么，问奈奎斯特吧。
还有就是output[0]对应着直流分量
5 _. R% p- n( o4 y& d6 ]
我只用了前几个意思意思。
& c& G8 D  p7 r3 G

1. cr4_fft_64_stm32(output,input,NPT);
   
2.         dsp_asm_powerMag(output);               
   
3. 
   
4.         LCD_ShowxNum(100,0,    OUTPUT_MAG[0]  ,7,16,0);        
   , k9 j9 e1 o! T: \+ S. F
5.         LCD_ShowxNum(100,16,   OUTPUT_MAG[1]  ,7,16,0);               
   
6.         LCD_ShowxNum(100,32,   OUTPUT_MAG[2]  ,7,16,0);               
   2 ~& C" ]. S- Q! n: b8 I' I4 W$ V  @
7.         LCD_ShowxNum(100,48,   OUTPUT_MAG[3]  ,7,16,0);               
   
8.         LCD_ShowxNum(100,64,   OUTPUT_MAG[4]  ,7,16,0);               
   * F; R: A; o* C3 G, {
9.         LCD_ShowxNum(100,80,   OUTPUT_MAG[5]  ,7,16,0);        
   
10.     THD=100*(float)sqrt((Voltage[4]*Voltage[4]+Voltage[2]*Voltage[2]+Voltage[3]*Voltage[3]+Voltage[5]*Voltage[5]))/(float)(Voltage[1]);
    

复制代码

————————————————
! U- s) w" y3 \' e; w6 d版权声明：ddydy111
$ ]7 ]) z  |! y7 w* e


2 }1 T' L* F- L; C1 L