寫在前面
20電賽整體感覺難度比之前小，本次程式設計上也沒有太多的難點，功能指標全部完成，程式實作了測量每種失真的情況下的THD的近似值，并且進行了程式拓展，實作了全自動的測量，以及顯示測量波形的波形圖，頻譜圖，根據題目要求，我們可以看出這次程式設計要用到FFT演算法，
我們的程式設計有兩個版本，一個版本是通過定時器進行采樣得到特定采樣率下的資料并保存在陣列里，然后進行傅里葉變換，另外一種就是通過定時器產生PWM波生成ADC的采樣時鐘，直接通過DMA保存資料然后進行傅里葉變換，
在理論計算下，
所以本文主要介紹我們實作的FFT的功能測驗驗證的程式，
題目
比賽指標要求
為什么需要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騙積分了，如果有點識訓希望能給個關注和點贊，
工程鏈接
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/175171.html
標籤：其他