2021-11-05 08:51:13 索煒達電子 1155
項目編號：E2145
文件大?。?8M
源碼說明：帶中文注釋
開發(fā)環(huán)境：C編譯器
簡要概述：
寫在前面
20電賽整體感覺難度比之前小，本次程序設計上也沒有太多的難點。功能指標全部完成，程序實現(xiàn)了測量每種失真的情況下的THD的近似值。并且進行了程序拓展，實現(xiàn)了全自動的測量，以及顯示測量波形的波形圖，頻譜圖。根據(jù)題目要求，我們可以看出這次程序設計要用到FFT算法。
我們的程序設計有兩個版本，一個版本是通過定時器進行采樣得到特定采樣率下的數(shù)據(jù)并保存在數(shù)組里，然后進行傅里葉變換，另外一種就是通過定時器產(chǎn)生PWM波生成ADC的采樣時鐘，直接通過DMA保存數(shù)據(jù)然后進行傅里葉變換。
所以本文主要介紹我們實現(xiàn)的FFT的功能測試驗證的程序。
比賽指標要求
為什么需要FFT?
任何連續(xù)測量的時域信號都可以表示為不同頻率的正弦波信號的無限疊加。以累加的方式來計算該信號中不同信號的頻率、振幅和相位。所以本次測量就必須要使用FFT算法。
知識科普：THD
總諧波失真表明功放工作時，由于電路不可避免的振蕩或其他諧振產(chǎn)生的二次，三次諧波與實際輸入信號疊加，在輸出端輸出的信號就不單純是與輸入信號完全相同的成分，而是包括了諧波成分的信號，這些多余出來的諧波成分與實際輸入信號的對比，用百分比來表示就稱為總諧波失真。一般來說，總諧波失真在500赫茲附近最小，所以大部分功放表明總諧波失真是用500赫茲信號做測試，但有些更嚴格的廠家也提供20－20000赫茲范圍內的總諧波失真數(shù)據(jù)?？傊C波失真在1%以下，一般耳朵分辨不出來，超過10%就可以明顯聽出失真的成分。這個總諧波失真的數(shù)值越小，音色就更加純凈。一般產(chǎn)品的總諧波失真都小于1%@500Hz，但這個數(shù)值越小，表明產(chǎn)品的品質越高。
所以在進行測試前我們就要先有個概念
對于信號源輸出的1k的正弦信號，總諧波失真的近似值越小，表示程序更精準，基本在1.0%以內。
對于信號源輸出的1k的方波信號，總諧波失真的近似值大約是0.3887（前5次諧波計算的近似值）
失真度測試儀測量的結果：
正弦波
方波
這里解釋下為啥方波測量出來的是44.26%，這里先給出方波的傅里葉變換式子
因為對于近似值來說方波取前五次傅里葉變換的值就是0.3887
計算到前7次時候
MTLAB仿真測試
所以根據(jù)已有的知識，進行下MATLAB仿真測試
clf;fs=10240; %采樣頻率
Ndata=1024; %數(shù)據(jù)長度
N=1024; %FFT的數(shù)據(jù)長度
n=0:Ndata-1;
t=n/fs; %數(shù)據(jù)對應的時間序列
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弧度
//復數(shù)的交換
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;
}
}
//求所有復數(shù)的模
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);
}
}
//求復數(shù)的和
void c_plus(complex a,complex b,complex *c)
{
c->real = a.real + b.real;
c->imag = a.imag + b.imag;
}
//求復數(shù)的差
void c_sub(complex a,complex b,complex *c)
{
c->real = a.real - b.real;
c->imag = a.imag - b.imag;
}
//求復數(shù)的積
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;
}
//求復數(shù)的商
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;
/*----計算分解的級數(shù)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級每個分組所含節(jié)點數(shù)
lb=la/2; //lb代表第m級每個分組所含碟形單元數(shù)
//同時它也表示每個碟形單元上下節(jié)點之間的距離
/*----碟形運算----*/
for(l=1;l<=lb;l++)
{
r=(l-1)*pow(2,M-m);
for(n=l-1;n<N-1;n=n+la) //遍歷每個分組，分組總數(shù)為N/la
{
lc=n+lb; //n,lc分別代表一個碟形單元的上、下節(jié)點編號
Wn_i(N,r,&wn,1);//wn=Wnr
c_mul(f[lc],wn,&t);//t = f[lc] * wn復數(shù)運算
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; /*每一級碟形運算間隔的點數(shù)*/
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 //復數(shù)類型
{
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);//復數(shù)加
void c_mul(complex a,complex b,complex *c) ;//復數(shù)乘
void c_sub(complex a,complex b,complex *c); //復數(shù)減法
void c_div(complex a,complex b,complex *c); //復數(shù)除法
void fft(int N,complex f[]);//傅立葉變換 輸出也存在數(shù)組f中
void ifft(int N,complex f[]); // 傅里葉逆變換
void c_abs(complex f[],float out[],int n);//復數(shù)數(shù)組取模
void FFT(struct compx *xin,int N);
#endif
main.c
程序是根據(jù)網(wǎng)上的程序更改參考的，用的是別人自己寫的FFT，把兩個人寫的放到了一起，大家可以根據(jù)需要自己選擇如何調用，后面會更新使用官方庫版本的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 //采樣點數(shù)
#define Fs 10240 //采樣頻率
#define F 10 //分辨率
/* Private variables ---------------------------------------------------------*/
/* Private function prototypes -----------------------------------------------*/
/* Private functions ---------------------------------------------------------*/
//FFT測試數(shù)據(jù)集 輸入數(shù)組
complex FFT_256PointIn[N];
//FFT測試數(shù)據(jù)集 輸出數(shù)組
float FFT_256PointOut[N/2];
//填入數(shù)組
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;
}
}
/******************************************************************
函數(shù)名稱:GetPowerMag()
函數(shù)功能:計算各次諧波幅值
參數(shù)說明:
備 注:先將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();//系統(tǒng)時鐘初始化
USART_Configuration();//串口1初始化
printf("這是一個FFT 測試實驗\r\n");
InitBufInArray();
fft(N,FFT_256PointIn);
printf("點數(shù) 頻率 幅值 實部 虛部\n");
//GetPowerMag();
dsp_g2_test();
while(1)
{
//檢查接收數(shù)據(jù)完成標志位是否置位
if(USART_GetFlagStatus(USART1, USART_IT_RXNE) != RESET)
{
//將接收到的數(shù)據(jù)發(fā)送出去，對USART_DR的讀操作可以將USART_IT_RXNE清零。
printf("%c",USART_ReceiveData(USART1));
}
}
}
/*********************************END OF FILE**********************************/
串口的截圖結果是正確的。
目錄│文件列表:
└ question-2020-e
│ demo.m
│ DMA fft（畫圖 自動 繼電器控制 全部顯示最終版本）.zip
│ DMA fft（畫圖 自動 繼電器控制-0.3版本）.zip
│ DMA fft（畫圖 自動-0.2版本）.zip
│ fft分度值10-定時器版本（頻譜圖 波形圖 THD）.zip
│ STM32的FFT官方庫.zip
│ 數(shù)據(jù)測試.zip
└ 電路設計部分.7z
