本文开始前我们先简单的介绍一下什么是THD?

" ~3 e5 `7 W+ }7 u; h! C总谐波失真（Total Harmonic Distortion，简称THD）是衡量信号失真程度的重要指标，特别是在音频、电力电子设备和通信系统等领域。THD表示一个信号中谐波分量（除了基波之外的频率分量）相对于基波的比例。具体来说，THD是所有谐波分量的平方和与基波分量平方的比值的平方根，通常用百分比表示。
$ h1 }& M3 D- L
公式表示为：


0 M" X8 G5 u4 R3 s! n0 J5 M# P/ P其中的V1是基波幅度，V2,V3……Vn是谐波幅度。
$ w  L$ F' Q8 o0 D0 C
. _/ o9 a& E- V, m1 K这个公式很好理解，一个信号的谐波幅度平方和占谐波频率的百分比。
; H9 `, Z9 T3 G1 V
那么如何获取基波幅度和谐波幅度，就是我们的FFT需要进行的事情。

+ |6 L6 w$ x( f7 K: T2 E/ C利用FFT我们可以获取信号的频率谱。



例如一段信号方波信号经过FFT之后，这些突出的方波信号就是谐波的幅值，而第一段（不包括第一个点）的则是其基波幅度。但是这里不包含第一个点，因为那个点是直流的幅度（两倍）。

因此我们只需要进行FFT之后统计谐波的平方和即可获得我们的THD。

快速傅里叶变换
9 w; h) `9 I  Q+ G& _: p- G: c关于FFT我们在之前的文章中有过介绍，我们在Keil中可以利用DSP库来进行快速的FFT（手搓的FFT算法效率通常很低）。我们在之前的文章中介绍过如何使用ADC+DMA的快速采样配合DSP库进行FFT。


# e, G1 o2 z  O+ \: O+ n" r
STM32的DMA采样+FFT时域分析（STM32F407）
( E' O+ }" p  O4 n0 f# y7 s& W
) k, C2 q. r/ Q, l. I) o' L我们同样的在这篇文章的基础上进行。
+ V" w5 ~% T  d7 ^! W4 t+ n# l5 R  o

1. #define ADCLenth 1024
   ; {, J  k1 w* I, L
2. uint16_t ADCValue[ADCLenth];
   7 B/ S3 C, C3 V4 d  ~  m+ b& v$ ^
3. arm_cfft_radix4_instance_f32 scfft;//定义scfft结构体
   ( Y( \: ^$ U5 [# @' U0 p; P+ w
4. float FFT_InputBuf[ADCLenth*2];  //FFT输入数组
   
5. float FFT_OutputBuf[ADCLenth];  //FFT输出数组

复制代码

3 i3 A3 v$ n/ j# t我们定义一个FFT的结构体变量，以及存放ADC采样结果的数组还有FFT的输入输出数组，这里的输入数组之所以长度翻倍是因为对于FFT的输入来说是一个复数即包含实部和虚部，因此长度是两倍。
; g7 k/ q) b* x/ w/ k1 E0 H

1.         for(int i=0; i < ADCLenth; i++)
   
2.         {
   ) [' Q& G. W1 d
3.           FFT_InputBuf[2*i]=(ADCValue[i] )*3.3/4096; //实部
   
4.           FFT_InputBuf[2*i+1]=0;          //虚部
   7 A8 e7 t! |) O# R) v) y3 p+ ]" `
5.         }
   
6.         arm_cfft_radix4_f32(&scfft,FFT_InputBuf);  
   
7.         arm_cmplx_mag_f32(FFT_InputBuf,FFT_OutputBuf,ADCLenth);  //取模得幅值
   ) z0 \& I) y( N3 Z5 u- K- d
8.         
   * ^% j8 R( X' s! R
9. 
   
10.         for(int i = 2;i<ADCLenth/2;i++)
    ! o  ~+ B0 O9 v% f$ e+ Q1 e
11.         {
    
12.           if(FFT_OutputBuf[i]>maxValue)
    7 j6 |5 U3 |/ A3 l  T6 _3 R
13.           {
    6 {( ]) `% G  r& [
14.          
    0 H$ R# F& `* K% [& `
15.             maxValue = FFT_OutputBuf[i];
    : ?- N  R) z  z
16.             max = i;
    
17.          
    
18.           }
    - L% B+ P  }7 s$ `( D7 |
19.         }

复制代码

一轮ADC采样结束之后，我们将其的实部信号和虚部信号（0）存放FFT的输入数组，之后执行快速傅里叶变换获得FFT的模值，模值存放在FFT_OutputBuf中。
) c, A# l( e0 b" H$ |
+ C! Y+ D7 {! s- B. }- X; H我们通过一个比较循环来寻找FFT结果的最大值，这里我们忽略前几个元素尤其是索引为0的位置，因为他是直流分量的模值。
* L; c  S& C& L: l  S; J/ ^4 t- \


! N  l4 q9 j) K* ?/ C0 l6 G之后我们就需要统计各个谐波的幅度。
& o, S7 ?, _7 X7 {, _" M8 V

1. float FindMax(float * fft,int index,int wind)
     m/ H, q- v1 U' ^4 }7 Z" [% ?
2. {
   ! ]5 a$ [% i9 T$ ^$ z1 a4 t
3.   int max = 0;
   * A8 @& E1 Q8 i; @  D7 F, F" ~! ^5 W; X
4.   for(int i = index-wind ;i <index+wind;i++)
   6 A& t& ~9 t; E) ?
5.   {
   0 f' q3 _2 g. O5 }" ^" X+ V
6.     if(fft[i]>max)
   
7.     {
   
8.       max =fft[i];
   % r" r3 i+ `$ I" ]" N$ N0 K3 i
9.     }
   
10.    
    
11.   }
    
12.   return max;
    
13.   
    
14. }

复制代码

* {7 ?: N2 O% `: f: R我们定义一个函数来寻找某点附近的最大值，这里之所以要寻找某点附近的最大值在这里说明一下。
5 T4 g: W1 a! `0 N2 o) X
FFT的索引和频率有关系，每个索引对应的频率之差为：采样频率/采样长度。我们的采样长度是1024，而采样频率是20kHZ,因此索引差为对应的频率差为19.53HZ，因为根据计算，1000HZ的频率对应的索引为51.2而由于索引只能是整数，因此在频谱上基波的最大值并不会是51.2而是51，因为int类型会抛弃小数。

所以我们实际统计的是984.3HZ或者1015.56HZ的频率，这样子我们计算谐波幅度进行翻倍的时候并不会是刚刚好好的对应的2KHZ的频率，而是其最大幅度会发生便宜。

/ g! i) V; d# E/ J) t( n! U因此我们因为是寻找我们认为的谐波索引的附近最大值。

1.         for(int i = max;i<ADCLenth/2;i+=max)
   ( q# ?0 f, Q+ b  d. f
2.         {
   - n/ x8 y3 N$ Y' w+ X! j
3.           maxW = FindMax(FFT_OutputBuf,i,max*0.3);
   
4.           all =all +  maxW*maxW;
   , l7 C6 m2 j4 A. t' e( \, O0 D
5.         }
   7 g  `! @( D  v; `+ }" t
6.         all = all - maxValue*maxValue;
   
7.         
   / u* J! o; }9 v- n; l5 y: i
8.         DHT = sqrtf(all)/maxValue;
   8 u2 N* ?" k- b; [/ t4 C' }
9.         printf("DHT:%f%c\r\n",DHT*100,'%');
   3 ?' n. y! h% S5 K& Y
10.         all = 0;
    
11.         max = 2;
      f/ a0 A& O; P9 I5 M+ ]
12.         maxValue = 0;

复制代码

这样子就是统计我们的THD的值，之后将其打印，这里计算谐波幅度的时候我把基波的幅度也加了上去，我们将其去除。

4 Y7 d7 _, q; E
0 t( e. \2 P5 [+ ]$ S$ m5 a$ f
* D0 Y& R, z- V9 c. N) O可以看到外面的THD（图中打错了）计算外面的方波频率的THD在40%左右。
( s* o4 Q" T( }' _- e0 b8 f我们利用MatLab生成一下方波信号之后用FFT来看一下理论结果。

0 v' X& D0 b# i
; ^- `& n! u2 }- I8 i. A& F
5 @" G3 P) p, J* ?1 Z

3 k( m+ N3 D% d* _1 P5 }其统计到5次谐波对应的值计算出的THD的值43%，这和我们的计算结果也很接近了。
* A2 e: `" L7 Y; X( a" s+ _我们同样的看一下手机上的显示内容。
' N! C4 {! u: `# ~$ \" ]
9 H1 P$ g, A$ x7 O& S

5 l7 B" _7 l. d5 Z0 R+ Y9 L总谐波失真为42.7%
) `: N/ K$ A6 r' c. ~: \
% k1 c5 s8 f8 Q7 o' s$ p* P. T
转载自：电路小白
7 P6 ~8 Q2 |; K( {* s  i如有侵权请联系删除

, V2 }6 `0 t5 \( P2 `( S( V
8 Q$ i3 j, }' j. X
  N2 z: E9 i9 [* v% J