相关文章推荐
逆袭的领结  ·  深度学习_卡拉比丘流形的博客-CSDN博客·  1 年前    · 
直爽的遥控器  ·  java - Is it required ...·  2 年前    · 
任性的柚子  ·  git提交忽略指定文件_yzh_134698 ...·  2 年前    · 
重感情的金鱼  ·  用户对问题“如何使用copy_from将带日 ...·  2 年前    · 
Code  ›  2020电赛E题--非线性失真器程序设计--01--算法仿真与STM32FFT数据验证（附工程源码+gitee链接）-阿里云开发者社区
非线性 电路仿真 fft 谐波失真
http://developer.aliyun.com:443/article/1053009
要出家的煎饼果子
2 年前
本文内容由阿里云实名注册用户自发贡献，版权归原作者所有，阿里云开发者社区不拥有其著作权，亦不承担相应法律责任。具体规则请查看《 阿里云开发者社区用户服务协议 》和 《 阿里云开发者社区知识产权保护指引 》。如果您发现本社区中有涉嫌抄袭的内容，填写 侵权投诉表单 进行举报，一经查实，本社区将立刻删除涉嫌侵权内容。

写在前面

20电赛整体感觉难度比之前小，本次程序设计上也没有太多的难点。功能指标全部完成，程序实现了测量每种失真的情况下的THD的近似值。并且进行了程序拓展，实现了全自动的测量，以及显示测量波形的波形图， 频谱图 。根据题目要求，我们可以看出这次程序设 计要用到FFT算法。

我们的程序设计有两个版本，一个版本是通过定时器进行 采样 得到特定采样率下的数据并保存在数组里，然后进行傅里叶变换，另外一种就是通过定时器产生PWM波生成ADC的采样时钟，直接通过DMA保存数据然后进行傅里叶变换。

在理论计算下，

所以本文主要介绍我们实现的 FFT 的功能测试验证的程序。

题目

image.png

比赛指标要求

image.png

为什么需要FFT?

任何连续测量的时域信号都可以表示为不同频率的正弦波信号的无限叠加。以累加的方式来计算该信号中不同信号的频率、振幅和相位。所以本次测量就必须要使用FFT算法。

原理部分就引用下别人的帖子大家自行查看：

超详细易懂FFT（快速傅里叶变换）及代码实现

知识科普：THD

总谐波失真表明功放工作时，由于电路不可避免的振荡或其他谐振产生的二次，三次谐波与实际输入信号叠加，在输出端输出的信号就不单纯是与输入信号完全相同的成分，而是包括了谐波成分的信号，这些多余出来的谐波成分与实际输入信号的对比，用百分比来表示就称为总谐波失真。一般来说，总谐波失真在500赫兹附近最小，所以大部分功放表明总谐波失真是用500赫兹信号做测试，但有些更严格的厂家也提供20－20000赫兹范围内的总谐波失真数据。总谐波失真在1%以下，一般耳朵分辨不出来，超过10%就可以明显听出失真的成分。这个总谐波失真的数值越小，音色就更加纯净。一般产品的总谐波失真都小于1%@500Hz，但这个数值越小，表明产品的品质越高。

所以在进行测试前我们就要先有个概念
对于信号源输出的1k的正弦信号，总谐波失真的近似值越小，表示程序更精准，基本在1.0%以内。
对于信号源输出的1k的方波信号，总谐波失真的近似值大约是0.3887（前5次谐波计算的近似值）

失真度测试仪测量的结果：

正弦波

image.png

方波

image.png

这里解释下为啥方波测量出来的是44.26%，这里先给出方波的 傅里叶 变换式子

image.png

因为对于近似值来说方波取前五次傅里叶变换的值就是0.3887

image.png

计算到前7次时候

image.png

MTLAB仿真测试

所以根据已有的知识，进行下MATLAB仿真测试 

clf;fs=10240; %采样频率
Ndata=1024; %数据长度
N=1024; %FFT的数据长度
n=0:Ndata-1;
t=n/fs;   %数据对应的时间序列
x=0.5*sin(2*pi*1000*t)+1;   %时间域信号
%subplot(2,2,4),plot(t,x);
subplot(2,2,2),plot(t,x,'.--');
y=fft(x,N);   %信号的Fourier变换
mag=abs(y);    %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');





    
title('Ndata=10240 Nfft=1024');grid on;

image.png

这里仿真显示的频谱图和我们的代码模拟给出的输入信号是相同的所以大致可以按照这个傅里叶变换的标准进行编写代码。之所以这里画出采样的波形图是因为后面我们程序要画波形图，所以这里就测试了下。理想波形的总谐波失真计算没有意义所以就不进行计算。

STM32测试程序：

FFT.C

这里的FFT也是找到的别人写好的程序，所以就不做详细讲解了（能力有限） 

#include "math.h"
#include "fft.h"
//精度0.0001弧度
//复数的交换 
void conjugate_complex(int n,complex in[],complex out[])
  int i = 0;
  for(i=0;i<n;i++)
    out[i].imag = -in[i].imag;
    out[i].real = in[i].real;
//求所有复数的模
void c_abs(complex f[],float out[],int n)
  int i = 0;
  float t;
  for(i=0;i<n;i++)
    t = f[i].real * f[i].real + f[i].imag * f[i].imag;
    out[i] = sqrt(t);
//求复数的和 
void c_plus(complex a,complex b,complex *c)
  c->real = a.real + b.real;
  c->imag = a.imag + b.imag;
//求复数的差  
void c_sub(complex a,complex b,complex *c)
  c->real = a.real - b.real;
  c->imag = a.imag - b.imag;    
//求复数的积
void c_mul(complex a,complex b,complex *c)





    
  c->real = a.real * b.real - a.imag * b.imag;
  c->imag = a.real * b.imag + a.imag * b.real;  
//求复数的商 
void c_div(complex a,complex b,complex *c)
  c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);
  c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);
#define SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr
void Wn_i(int n,int i,complex *Wn,char flag)
  Wn->real = cos(2*PI*i/n);
  if(flag == 1)
  Wn->imag = -sin(2*PI*i/n);
  else if(flag == 0)
  Wn->imag = -sin(2*PI*i/n);
//傅里叶变化
void fft(int N,complex f[])
  complex t,wn;//中间变量
  int i,j,k,m,n,l,r,M;
  int la,lb,lc;
  /*----计算分解的级数M=log2(N)----*/
  for(i=N,M=1;(i=i/2)!=1;M++); 
  /*----按照倒位序重新排列原信号----*/
  for(i=1,j=N/2;i<=N-2;i++)
    if(i<j)
      t=f[j];
      f[j]=f[i];
      f[i]=t;
    k=N/2;
    while(k<=j)
      j=j-k;
      k=k/2;





    
    j=j+k;
  /*----FFT算法----*/
  for(m=1;m<=M;m++)
    la=pow(2,m); //la=2^m代表第m级每个分组所含节点数     
    lb=la/2;    //lb代表第m级每个分组所含碟形单元数
                 //同时它也表示每个碟形单元上下节点之间的距离
    /*----碟形运算----*/
    for(l=1;l<=lb;l++)
      r=(l-1)*pow(2,M-m);   
      for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la
        lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号     
        Wn_i(N,r,&wn,1);//wn=Wnr
        c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算
        c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr
        c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr
//傅里叶逆变换
void ifft(int N,complex f[])
  int i=0;
  conjugate_complex(N,f,f);
  fft(N,f);
  conjugate_complex(N,f,f);
  for(i=0;i<N;i++)
    f[i].imag = (f[i].imag)/N;
    f[i].real = (f[i].real)/N;
struct compx EE(struct compx b1,struct compx b2)
    struct compx b3;
    b3.real = b1.real*b2.real-b1.imag*b2.imag;
    b3.imag = b1.real*b2.imag+b1.imag*b2.real;
    return (b3);





    
void FFT(struct compx *xin,int N)
int f,m,LH,nm,i,k,j,L;
double p , ps ;
int le,B,ip;
float pi;
struct compx w,t;
LH=N/2; 
for(m=1;(f=f/2)!=1;m++){;}  /*2^m=N*/
for(L=m;L>=1;L--)    /*这里和时域的也有差别*/
le=pow(2,L);
B=le/2; /*每一级碟形运算间隔的点数*/
pi=3.14159;
 for(j=0;j<=B-1;j++)
   p=pow(2,m-L)*j;
   ps=2*pi/N*p;
   w.real=cos(ps);
   w.imag=-sin(ps);
   for(i=j;i<=N-1;i=i+le)
      ip=i+B;  
      t=xin[i];
      xin[i].real=xin[i].real+xin[ip].real;
      xin[i].imag=xin[i].imag+xin[ip].imag;  
      xin[ip].real=xin[ip].real-t.real;
      xin[ip].imag=xin[ip].imag-t.imag;     
      xin[ip]=EE(xin[ip],w);
/*变址运算*/
nm=N-2;   
j=N/2;
for(i=1;i<=nm;i++)
if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=LH;





    
while(j>=k){j=j-k;k=k/2;}
j=j+k;
}

FFT.H 

#ifndef __FFT_H__
#define __FFT_H__
typedef struct complex //复数类型
  float real;       //实部
  float imag;       //虚部
}complex;
struct compx
    double real;
    double imag;
#define PI 3.1415926535897932384626433832795028841971
void conjugate_complex(int n,complex in[],complex out[]);
void c_plus(complex a,complex b,complex *c);//复数加
void c_mul(complex a,complex b,complex *c) ;//复数乘
void c_sub(complex a,complex b,complex *c); //复数减法
void c_div(complex a,complex b,complex *c); //复数除法
void fft(int N,complex f[]);//傅立叶变换 输出也存在数组f中
void ifft(int N,complex f[]); // 傅里叶逆变换
void c_abs(complex f[],float out[],int n);//复数数组取模
 void FFT(struct compx *xin,int N);
#endif

main.c

程序是根据网上的程序更改参考的，用的是别人自己写的FFT，把两个人写的放到了一起，大家可以根据需要自己选择如何调用，后面会更新使用官方库版本的FFT的代码版本。串口部分就使用串口1就行，如果是正点的板子程序改写是默认打开了串口1的。 

/* Includes ------------------------------------------------------------------*/
#include "main.h"





    
#include "usart.h"
#include "fft.h"
#include <math.h>
/* Private typedef -----------------------------------------------------------*/
/* Private define ------------------------------------------------------------*/
/* Private macro -------------------------------------------------------------*/
#define  N    1024          //采样点数
#define  Fs   10240        //采样频率
#define  F    10          //分辨率
/* Private variables ---------------------------------------------------------*/
/* Private function prototypes -----------------------------------------------*/
/* Private functions ---------------------------------------------------------*/
//FFT测试数据集 输入数组
complex  FFT_256PointIn[N];
//FFT测试数据集 输出数组
float   FFT_256PointOut[N/2];                                       
//填入数组          
double result[N];
struct compx s[N];
void InitBufInArray()
 unsigned short i;
 for(i=0; i<N; i++)    
       FFT_256PointIn[i].real  = 1 * sin(2*PI * i * 1000.0 / Fs) 
           FFT_256PointIn[i].imag = 0;
/******************************************************************
函数名称:GetPowerMag()
函数功能:计算各次谐波幅值
参数说明:
备　　注:先将FFT_256PointIn分解成实部(X)和虚部(Y)，





    
         然后计算幅值:(sqrt(X*X+Y*Y)*2/N
         然后计算相位:atan2(Y/X)
*******************************************************************/
void GetPowerMag()
    unsigned short i;
      float  X,Y,P,Mag;
        c_abs(FFT_256PointIn,FFT_256PointOut,N/2);
    for(i=0; i<N/2; i++)
              X = FFT_256PointIn[i].real/N;    //计算实部
              Y = FFT_256PointIn[i].imag/N;    //计算虚部
              Mag = FFT_256PointOut[i]*2/N;    //计算幅值
              P = atan2(Y,X)*180/PI;           //计算相位
                printf("%d      ",i);
                printf("%d      ",F*i); 
                printf("%f      \r\n",Mag);
//            printf("%f      ",P);
//              printf("%f      ",X);
//              printf("%f      \r\n",Y);           
void dsp_g2_test()
  u16 i=0;
  for(i=0;i<N;i++)
    s[i].real = 32000 * sin(PI*2*i*(50.0f/Fs));
    s[i].real+= 16000 * sin(PI*2*i*(550.0f/Fs));
    s[i].real+= 9000 * sin(PI*2*i*(1150.0f/Fs));
        s[i].real+= 6000 * sin(PI*2*i*(2100.0f/Fs));
        s[i].real+= 4000 * sin(PI*2*i*(5000.0f/Fs));
    s[i].imag=0;
  FFT(s,N);
  for(i=0;i<N/2;i++)
        if(i==0)





    
            result[i] = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag)/N;
            result[i] = sqrt(s[i].real * s[i].real + s[i].imag * s[i].imag)*2/N;
        printf("%d      ",i);
        printf("%d      ",10*i);
        printf("%f      \r\n",result[i]);
    //if(result[i] > 10)
    //printf("%4d,%4d,%ld\n",i,(u16)((double)i*Fs/NPT),(u32)result[i]);
  * @brief  串口打印输出
  * @param  None
  * @retval None
int main(void)
    int i,t;
    SystemInit();//系统时钟初始化
    USART_Configuration();//串口1初始化
    printf("这是一个FFT 测试实验\r\n");  
  InitBufInArray(); 
    fft(N,FFT_256PointIn);
    printf("点数   频率  幅值   实部  虚部\n"); 
    //GetPowerMag();
    dsp_g2_test();
    while(1)
        //检查接收数据完成标志位是否置位   
        if(USART_GetFlagStatus(USART1, USART_IT_RXNE) != RESET)
        //将接收到的数据发送出去，对USART_DR的读操作可以将USART_IT_RXNE清零。
        printf("%c",USART_ReceiveData(USART1));
/*********************************END OF FILE**********************************/

串口的截图结果是正确的。

image.png

免费开源，大家参考，不通过csdn骗积分了，如果有点收获希望能给个关注和点赞。

工程链接

image.png

 
推荐文章
逆袭的领结  ·  深度学习_卡拉比丘流形的博客-CSDN博客
1 年前
直爽的遥控器  ·  java - Is it required to use orElse(null) in Streams? - Stack Overflow
2 年前
任性的柚子  ·  git提交忽略指定文件_yzh_1346983557的博客-CSDN博客
2 年前
重感情的金鱼  ·  用户对问题“如何使用copy_from将带日期时间的csv加载到postgresql中”的回答 - 问答 - 腾讯云开发者社区-腾讯云
2 年前
今天看啥   ·   Py中国   ·   codingpro   ·   小百科   ·   link之家   ·   卧龙AI搜索
删除内容请联系邮箱 2879853325@qq.com
Code - 代码工具平台
© 2024 ~ 沪ICP备11025650号