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