写在前面
20电赛整体感觉难度比之前小，本次程序设计上也没有太多的难点。功能指标全部完成，程序实现了测量每种失真的情况下的THD的近似值。并且进行了程序拓展，实现了全自动的测量，以及显示测量波形的波形图，频谱图。根据题目要求，我们可以看出这次程序设计要用到FFT算法。
我们的程序设计有两个版本，一个版本是通过定时器进行采样得到特定采样率下的数据并保存在数组里，然后进行傅里叶变换，另外一种就是通过定时器产生PWM波生成ADC的采样时钟，直接通过DMA保存数据然后进行傅里叶变换。
在理论计算下，
所以本文主要介绍我们实现的FFT的功能测试验证的程序。
题目
比赛指标要求
为什么需要FFT?
任何连续测量的时域信号都可以表示为不同频率的正弦波信号的无限叠加。以累加的方式来计算该信号中不同信号的频率、振幅和相位。所以本次测量就必须要使用FFT算法。
原理部分就引用下别人的帖子大家自行查看：
知识科普：THD
总谐波失真表明功放工作时，由于电路不可避免的振荡或其他谐振产生的二次，三次谐波与实际输入信号叠加，在输出端输出的信号就不单纯是与输入信号完全相同的成分，而是包括了谐波成分的信号，这些多余出来的谐波成分与实际输入信号的对比，用百分比来表示就称为总谐波失真。一般来说，总谐波失真在500赫兹附近最小，所以大部分功放表明总谐波失真是用500赫兹信号做测试，但有些更严格的厂家也提供20－20000赫兹范围内的总谐波失真数据。总谐波失真在1%以下，一般耳朵分辨不出来，超过10%就可以明显听出失真的成分。这个总谐波失真的数值越小，音色就更加纯净。一般产品的总谐波失真都小于1%@500Hz，但这个数值越小，表明产品的品质越高。
所以在进行测试前我们就要先有个概念
对于信号源输出的1k的正弦信号，总谐波失真的近似值越小，表示程序更精准，基本在1.0%以内。
对于信号源输出的1k的方波信号，总谐波失真的近似值大约是0.3887（前5次谐波计算的近似值）
失真度测试仪测量的结果：
正弦波
方波
这里解释下为啥方波测量出来的是44.26%，这里先给出方波的傅里叶变换式子
因为对于近似值来说方波取前五次傅里叶变换的值就是0.3887
计算到前7次时候
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;
这里仿真显示的频谱图和我们的代码模拟给出的输入信号是相同的所以大致可以按照这个傅里叶变换的标准进行编写代码。之所以这里画出采样的波形图是因为后面我们程序要画波形图，所以这里就测试了下。理想波形的总谐波失真计算没有意义所以就不进行计算。
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; f=N; 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) +1; 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; else 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**********************************/
串口的截图结果是正确的。
免费开源，大家参考，不通过csdn骗积分了，如果有点收获希望能给个关注和点赞。
