目录
- 计算公式
- 波形生成
- DAC波形频率
- ADC采样时间
- 离散傅里叶变换DFT
- FFT
计算公式
DAC、ADC、FFT之间有些参数环环相扣,所以先整合一下公式。
1.系统时钟周期72MHZ。
2.定时器的单个时钟周期。
3.定时器的触发周期。
4.正弦波一个周期的时间,其中N为正弦波一个周期的点数。
5.正弦波的频率。
6.adc的采样频率,其中n为一个波形周期的采样点,fsin为被采样的波形的频率。
7.为设置adc的采样频率,要根据这个公式去配置触发adc的定时器。
8.进行傅里叶变换后,fft输出数组下标对应的频率。其中i为数组下标,fadc为adc的采样频率,fftnum为fft计算的点数。
波形生成
x取值范围[0,2π]
y=sin(x)取值范围[-1,1]
y=six(x)+1取值范围[0,2]
DAC输出电压范围[0,Vmax]
将y取值范围扩大到DAC输出电压范围,只需y=((six(x)+1)/2)*Vmax
周期2π,波形一个周期的点数为N,两点间距2π/N
下面这个代码,增加了一个DAClength为的是与DAC的DMA=normal配合,使DAC输出几个周期波形后,停止输出,满足一些特定电路的需求(有的电路起始阶段需要几个波形之后就停止在一个高电位,并不需持续输出波形)。
如果仅生成一个波形周期,那么可以无视DAClength参数和下面的复制周期。
/*** 生成正弦波数据点函数* @param NPoints 一个周期内的点数* @param DAClength 目的输出点数总数,DAClength为NPoints整数倍* @param VMaxRange 输出的电压最大值,取值范围0~3.3V* @param SineWaveTable 存放生成的数据点*/
void SineWaveGen(uint32_t NPoints, uint32_t DAClength, float VMaxRange, uint16_t* SineWaveTable)
{
#ifndef PI
#define PI 3.14159265358979323846
#endifint i = 0;int j = 0; int k = DAClength/NPoints; //增加波形后的周期数double radian = 0; // 弧度double setup = 0; // 弧度和弧度之间的大小double voltage = 0; // 输出电压setup = (2 * PI) / NPoints; // 两点之间的间距while (i < NPoints){voltage = VMaxRange / 2.0 * (sin(radian) + 1.0); // 计算电压//printf("%d %lf\r\n",i,voltage);SineWaveTable[i] = (uint16_t)(voltage * 4095 / 3.3); // 电压转为DAC数值//printf("%d %d\r\n",i,SineWaveTable[i]);radian += setup; // 下一个点的弧度i++;}for(j=1; j<k;++j)//复制k-1个周期{for(i=0;i<NPoints;++i){SineWaveTable[NPoints*j+i]=SineWaveTable[i];//printf("%d %d\r\n",NPoints*j+i,SineWaveTable[NPoints*j+i]);}}}
DAC波形频率
设置波形一个周期的点数,会影响DAC输出波形频率。
如果波形一个周期128个点,10k频率。不改变定时器设置的话,波形一个周期256个点,输出波形频率就变成了5k。
如果按下面配置定时器触发DAC,触发频率为72M/141=510638HZ。
波形一个周期点数为128,那么波形频率为510638HZ/128=3,989HZ,约为4k。
如果编写如下代码,波形一个周期点数为128,但是复制了九个周期的波形到dac数组里。
并且设置DMA的数据传输数量(0至65535)为128*9,最终得到的波形仍然是4k。
如果配置DMA Mode为Normal的话,那么可以发现触发一次DAC,输出了九个周期的频率为4k的波形。
#define POINTS 128
#define DAC_length 1152
uint16_t SineWaveTable[DAC_length];
SineWaveGen(POINTS,DAC_length, 2, SineWaveTable);//points扩展到daclength
HAL_DAC_Start_DMA(&hdac, DAC_CHANNEL_1, (uint32_t *)SineWaveTable, DAC_length, DAC_ALIGN_12B_R);
ADC采样时间
ADC使用若干个ADC_CLK周期对输入电压采样,采样周期数目可以通过ADC_SMPR1和ADC_SMPR2寄存器中的SMP[2:0]位更改。每个通道可以分别用不同的时间采样 。
总转换时间 :TCONV = 采样时间+ 12.5个周期
设置ADC输入时钟为12Mhz,那么1个ADC周期占用的时间=1 / 12MHZ = 0.0833334 uS
如果设置采样时间为1.5个周期,那么一次采样总的时间 = 采样时间 + 12.5个周期 = 1.5周期 + 12.5周期 = 14周期 = 14 * 0.0833334 = 1.166667 uS
两次采样间隔时间 = 1.166667 uS (ADCCLK为12MHZ时的最小采样间隔时间)
下面是ADC转换的时序图,可知,定时器触发adc转换的时间间隔需要大于adc两次采样间隔时间。
如果按照下面这样配置触发adc的定时器,72M/141=510638HZ
转换成时间就是1.95us,大于adc两次采样间隔时间,所以ok。
离散傅里叶变换DFT
一开始其实是打算手写一个DFT,但是实际用的时候无法满足单片机性能需求。单片机里面的算法,最好只有加减乘,不能有除。
离散傅里叶变换:
N为时域离散信号的点数,n为时域离散信号的编号(取值范围为0 ~ N-1),m为频域信号的编号(取值范围为0 ~ N-1),频域信号的点数也为N。
离散傅里叶变换的输入为N个离散的点(时域信号),输出为N个离散的点(频域信号,频域信号的每个点都用一个复数表示)。
那么可以根据以上公式写一个dft:
#include<iostream>
#include<cmath>
using namespace std;
double a[128] = {2.028,
2.038,
2.041,
2.032,
2.020,
2.026,
2.070,
2.070,
2.060,
2.061,
2.042,
2.037,
2.072,
2.080,
2.064,
2.063,
2.037,
2.032,
2.060,
2.051,
2.038,
2.034,
2.003,
1.990,
2.007,
2.005,
1.988,
1.978,
1.942,
1.936,
1.951,
1.940,
1.914,
1.902,
1.862,
1.850,
1.858,
1.852,
1.822,
1.810,
1.760,
1.755,
1.767,
1.751,
1.725,
1.716,
1.659,
1.648,
1.659,
1.652,
1.625,
1.607,
1.560,
1.544,
1.566,
1.550,
1.525,
1.510,
1.469,
1.454,
1.481,
1.469,
1.440,
1.434,
1.389,
1.380,
1.406,
1.401,
1.383,
1.377,
1.337,
1.333,
1.365,
1.360,
1.345,
1.339,
1.309,
1.310,
1.345,
1.348,
1.342,
1.347,
1.316,
1.322,
1.359,
1.358,
1.352,
1.355,
1.338,
1.347,
1.393,
1.397,
1.395,
1.406,
1.388,
1.403,
1.452,
1.459,
1.464,
1.474,
1.460,
1.479,
1.522,
1.535,
1.547,
1.561,
1.545,
1.562,
1.611,
1.625,
1.633,
1.651,
1.640,
1.659,
1.707,
1.718,
1.728,
1.743,
1.727,
1.738,
1.792,
1.803,
1.798,
1.818,
1.807,
1.816,
1.865,
1.873,};
double b[128];
const double PI = acos(-1.0);
//定义一个结构体来描述一个复数
typedef struct {float re;// reallyfloat im;// imaginary
} complex,*pcomplex;//构建并初始化一个复数结构体
complex complexBuild(float re,float im) {complex cx;cx.re=re;cx.im=im;return cx;
}//复数加法
complex complexAdd(complex a,complex b) {complex ret;ret.re=a.re+b.re;ret.im=a.im+b.im;return ret;
}
//复数乘法
complex complexMult(complex a,complex b) {complex ret;ret.im=a.im*b.re+a.re*b.im; ret.re=a.re*b.re-a.im*b.im;return ret;
}void DFT(complex x[],complex X[],int N) {int k,n;complex Wnk;for (k=0; k<N; k++) {X[k].re=0;X[k].im=0;for (n=0; n<N; n++) {//带公式 Wnk.re=(float)cos(2*PI*k*n/N);Wnk.im=(float)-sin(2*PI*k*n/N);X[k]=complexAdd(X[k],complexMult(x[n],Wnk));}}
}int main() {complex samples[128],_out[128];double _out2[128];int i;for(int i=0; i<128; ++i) {b[i]=a[i];samples[i].re=b[i];samples[i].im=0;printf("%.3f\n",b[i]);//printf("%.3f\n",b[i]);}
// for (i=0; i<120; i++) {
// samples[i].re=i;
// samples[i].im=0;
// }DFT(samples,_out,128);//求幅值 for(i=0;i<128;++i){_out2[i]=sqrt(_out[i].re*_out[i].re+_out[i].im*_out[i].im);}
// for (i=0; i<120; i++) {
// if(i==0)
// printf("(%f,%f)\n",_out[i].re/128,_out[i].im);
// else
// printf("(%f,%f)\n",_out[i].re/64,_out[i].im);
// }//数据处理 for (i=0; i<128; i++) {if(i==0)printf("%f\n",_out2[i]/128);elseprintf("%f\n",_out2[i]/64);}
}
/*
int main() {//memset(b,0,sizeof(b));for(int i=0; i<120; ++i) {b[i]=a[i]/4096.0*3.3;printf("%.3f\n",b[i]);}}
*/
由dft变换后的幅度可以看出波形的一些特征。比如直流偏置为1.674602。
fi是,进行傅里叶变换后,fft输出数组下标对应的频率。其中i为数组下标,fadc为adc的采样频率,fftnum为fft计算的点数。
按照之前的设置,adc的采样频率为510638HZ,fft计算的点数为128,那么
i=1时,f1=510638HZ/128=3989HZ。约等于4k,刚好对应上dac输出的正弦波的频率。
FFT
最后选择使用dsp库里面的fft进行傅里叶变换。
下面是一部分核心代码。直接做1024个点的fft。
#include "arm_math.h"
#include "arm_const_structs.h"
#define FFT_LENGTH 1024
float fft_inputbuf[FFT_LENGTH * 2];
float fft_outputbuf[FFT_LENGTH];
//部分代码:for (int i = 0; i < FFT_LENGTH; i++){fft_inputbuf[i * 2] = adc1_buff[i] * 3.0 / 4095;//实部赋值,* 3 / 4096是为了将ADC采集到的值转换成实际电压//printf("%.4f\r\n",fft_inputbuf[i * 2]);fft_inputbuf[i * 2 + 1] = 0;//虚部赋值,固定为0.}arm_cfft_f32(&arm_cfft_sR_f32_len1024, fft_inputbuf, 0, 1);arm_cmplx_mag_f32(fft_inputbuf, fft_outputbuf, FFT_LENGTH); /*处理变换结果*/fft_outputbuf[0] /= FFT_LENGTH;for (int i = 1; i < FFT_LENGTH; i++)//输出各次谐波幅值{fft_outputbuf[i] /= FFT_LENGTH/2;}/*打印结果*/printf("FFT Result:\r\n");for (int i = 0; i < FFT_LENGTH; i++)//输出各次谐波幅值{printf("%d:\t%.4f\r\n", i, fft_outputbuf[i]);}
但是实际操作过程中,模拟电路设计的有问题,导致adc读取数据在一定范围向上偏斜。
所以改变思路,做128个点的fft,去掉最大最小取平均。
adc读取数据一直向上偏,是电路本身存在问题,如果用算法去抵消这个影响,其实并没有解决本质问题。
int bnum = FFT_LENGTH / POINTS;//倍数int inum = 1;//4k对应fft128下标float resultzhi[bnum];//直流float resultfen[bnum];//分流for(int j=0;j<bnum;++j){for(int i=0;i<POINTS;++i){fft_inputpoint[i*2] = adc1_buff[j*POINTS+i] * 3.0 / 4095;fft_inputpoint[i*2+1] = 0;}arm_cfft_f32(&arm_cfft_sR_f32_len128, fft_inputpoint, 0, 1);arm_cmplx_mag_f32(fft_inputpoint, fft_outputpoint, POINTS); resultzhi[j]=fft_outputpoint[0]/POINTS;resultfen[j]=fft_outputpoint[inum]/(POINTS/2);//fft_outputpoint[inum]/=POINTS/2;//resultfen[j]=fft_outputpoint[inum];}float zhi=0,fen=0;for(int i=1;i<bnum-1;++i)//一定是首位最值{zhi+=resultzhi[i];fen+=resultfen[i];}zhi=zhi/(bnum-2);fen=fen/(bnum-2);printf("%.4f\r\n",zhi);printf("%.4f\r\n",fen);
后面经过一系列测试,得到一系列非线性的公式,无论选用何种拟合手段,都无法满足精度需求,所以,只能继续修改模拟电路。
由此,感悟就是,测量得到两个值之间不是线性关系。此时很多人会从算法层面切入。
但是还有一种方式就是从硬件层面切入,使用合适的电路,让两个值之间是线性关系。