python 失真度测试 失真度测试仪的原理

将失真度仪和示波器放到一块是因为两个工程需要做的事情非常的相似，只有小部分问题不同。所以将两个放在一起说，两个的不同点我会标注出来。

信号系统中的“失真度”定义为全部谐波能量与基波能量之比的平方根值。失真有多种：谐波失真、互调失真、相位失真等等。我们通常所说的失真度的技术术语为总谐波失真，简称THD。

失真度的测量方法

1.基波抑制法

2.谐波分析

3.模拟法

4.数字化方法

数字化方法是指先通过将信号数字化并送入计算机，在由计算机计算出失真度的测量方法。根据失真度的计算方法可分为FFT法和曲线拟合法。

利用stm32肯定是采用FFT（快速傅里叶变换）进行测量。

一、THD基本定义公式：
THD的基本定义是功率比例关系，各次谐波功率值和与基波功率之比的百分比表示

二、THD等效关系公式：

三、THD近似测量公式：

在数字信号处理中，学习了FFT是时域一频域变换分析中最基本的方法之一，是非常重要的一个概念。这里对FFT就不再过多介绍。

失真度仪代码整体思路

利用STM32将一个正弦波采集64或者256或者1024个点，一定是要保证最少采满一个周期，之所以采64或者256或者1024点是因为，st官方提供了一个FFT变换的库，链接：https://pan.baidu.com/s/1rVO6uYRk7eL2yID2U8IF9Q 提取码：mhsr,现在官方已经下载不了了，我把这个库的文件分享给大家。

采好点后利用DMA传输（减轻CPU负担），做FFT变换，做完以后直接套用THD计算公式计算失真度。

示波器代码整体思路

示波器工作和失真度仪非常雷同，采点，采好后利用DMA传输进行FFT，然后得到一些想要的值，比如最大值，最小值，信号的频率等等。

需要注意的是，这里AD采集需要配置为定时器触发，这样的话可以定时器定时时间就是 adc的采样时间，采样频率>2*信号频率。后续再通过按键调整采样率，这样2021年电赛A题不就出来了（电赛用的TI板子）。

失真度仪：

adc配置
 #include "adc.h"
//初始化定时器                                                               
 void  Adc_Init(void)
 {     
     ADC_InitTypeDef ADC_InitStructure; 
     GPIO_InitTypeDef GPIO_InitStructure;    RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA |RCC_APB2Periph_ADC1, ENABLE );      //使能ADC1通道时钟
      RCC_ADCCLKConfig(RCC_PCLK2_Div6);   //设置ADC分频因子6 72M/6=12,ADC最大时间不能超过14M
    //PA6 作为模拟通道输入引脚                         
     GPIO_InitStructure.GPIO_Pin = GPIO_Pin_6;
     GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AIN;        //模拟输入
     GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;
     GPIO_Init(GPIOA, &GPIO_InitStructure);        ADC_DeInit(ADC1);  //复位ADC1,将外设 ADC1 的全部寄存器重设为缺省值
    ADC_InitStructure.ADC_Mode = ADC_Mode_Independent;    //ADC工作模式:ADC1工作在独立模式
     ADC_InitStructure.ADC_ScanConvMode = DISABLE;    //模数转换工作在单通道模式
     ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;    //模数转换工作在非连续转换模式
     ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_T2_CC2;    //转换由定时器2的通道2触发（只有在上升沿时可以触发）
     ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;    //ADC数据右对齐
     ADC_InitStructure.ADC_NbrOfChannel = 1;    //顺序进行规则转换的ADC通道的数目
     ADC_Init(ADC1, &ADC_InitStructure);    //根据ADC_InitStruct中指定的参数初始化外设ADCx的寄存器       ADC_Cmd(ADC1, ENABLE);    //使能指定的ADC1
     
     ADC_DMACmd(ADC1, ENABLE);    //ADC的DMA功能使能
     
     ADC_ResetCalibration(ADC1);    //使能复位校准  
      
     ADC_RegularChannelConfig(ADC1, ADC_Channel_6, 1, ADC_SampleTime_1Cycles5 );//ADC1通道6,采样时间为239.5周期     
      
     ADC_ResetCalibration(ADC1);//复位较准寄存器
      
     while(ADC_GetResetCalibrationStatus(ADC1));    //等待复位校准结束
     
     ADC_StartCalibration(ADC1);     //开启AD校准
  
     while(ADC_GetCalibrationStatus(ADC1));     //等待校准结束
  
     ADC_SoftwareStartConvCmd(ADC1, ENABLE);        //使能指定的ADC1的软件转换启动功能}

         

 ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_T2_CC2;    //转换由定时器2的通道2触发（只有在上升沿时可以触发）

一定要注意这句，配置为定时器触发。

将采集好的值作为一个数组调用以上函数进行FFT变换，然后套用公式计算失真度

//计算总谐振失真函数
 void GetTHD()
 {
     unsigned short i=20;
     float  Uo1,Uo2,Uo3,Uo4,Uo5;
     double THD,thd_fz=0,thd_fm=0;
     Uo1=Mag[100];
     Uo2=Mag[200];
     Uo3=Mag[300];
     Uo4=Mag[400];
     Uo5=Mag[500];
     thd_fm=Uo1;
     thd_fz=Uo2*Uo2 + Uo3*Uo3 + Uo4*Uo4 + Uo5*Uo5;
     thd_fz=sqrt(thd_fz);
     THD=thd_fz/thd_fm*100;
     sprintf((unsigned char *)temp1,"THD:%5lf%%",THD);
 }

接下来再说一下调整采样率，将按键设置为中断线，按下按键进入中断，在中断服务函数中，改变定时器定时时间就改变了采样率，在标准库里，有这样一个函数：

void TIM_PrescalerConfig(TIM_TypeDef* TIMx, uint16_t Prescaler, uint16_t TIM_PSCReloadMode)
 {
   /* Check the parameters */
   assert_param(IS_TIM_ALL_PERIPH(TIMx));
   assert_param(IS_TIM_PRESCALER_RELOAD(TIM_PSCReloadMode));
   /* Set the Prescaler value */
   TIMx->PSC = Prescaler;
   /* Set or reset the UG Bit */
   TIMx->EGR = TIM_PSCReloadMode;
 }

很显然，这个函数可以改变定时器的pre，即adc的采样率。

总结一下，adc使用定时器触发设置采样率，DMA传输进行FFT变换，可是ST官方提供的FFT库只有64点，256点和1024点。如果我想做32点的怎么办，如果我想做512点的怎么办，这时候官方的库就不在适用了，只能自己写一个FFT变换的函数，我写了一个蝶形变换的函数，可以自己设置点数，只要是2^n次方就行。

void Mfft(double pr[], double pi[], int n, int k, double fr[], double fi[])
//pr[]  输入实部
 //  pi[]  输入虚部  （时域时未进行傅里叶变换，没有虚部，此数组全部赋值为0）
 // n 采样点数
 //  k :  2的k次方=n      采样点数可以自定义
 //fr[]  变换实部
 //  fi[]  变换虚部{
     int it,m,is,i,j,nv,l0;
     double p,q,s,vr,vi,poddr,poddi;
     for (it=0; it<=n-1; it++)  //将pr[0]和pi[0]循环赋值给fr[]和fi[]
     {
         m=it;
         is=0;
         for(i=0; i<=k-1; i++)
         {
             j=m/2;
             is=2*is+(m-2*j);
             m=j;
         }
         fr[it]=pr[is];
         fi[it]=pi[is];
      }
      pr[0]=1.0;
      pi[0]=0.0;
      p=6.283185306/(1.0*n);
      pr[1]=cos(p); //w=e^-j2pi/n欧拉公式表
      pi[1]=-sin(p);     for (i=2; i<=n-1; i++)  //计算pr[]
      {
          p=pr[i-1]*pr[1];
          q=pi[i-1]*pi[1];
          s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
          pr[i]=p-q; pi[i]=s-p-q;
      }
      for (it=0; it<=n-2; it=it+2)
      {
           vr=fr[it];
           vi=fi[it];
           fr[it]=vr+fr[it+1];
           fi[it]=vi+fi[it+1];
           fr[it+1]=vr-fr[it+1];
           fi[it+1]=vi-fi[it+1];
      }
      m=n/2;
      nv=2;
      for (l0=k-2; l0>=0; l0--) //蝶形计算
      {
           m=m/2;
           nv=2*nv;
           for (it=0; it<=(m-1)*nv; it=it+nv)
               for (j=0; j<=(nv/2)-1; j++)
               {
                     p=pr[m*j]*fr[it+j+nv/2];
                     q=pi[m*j]*fi[it+j+nv/2];
                     s=pr[m*j]+pi[m*j];
                     s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
                     poddr=p-q;
                     poddi=s-p-q;
                     fr[it+j+nv/2]=fr[it+j]-poddr;
                     fi[it+j+nv/2]=fi[it+j]-poddi;
                     fr[it+j]=fr[it+j]+poddr;
                     fi[it+j]=fi[it+j]+poddi;
                }
        }
        for (i=0; i<=n-1; i++)
        {
            pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);  //计算幅值
        }       return;
 }