DAC、ADC、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);

后面经过一系列测试，得到一系列非线性的公式，无论选用何种拟合手段，都无法满足精度需求，所以，只能继续修改模拟电路。

由此，感悟就是，测量得到两个值之间不是线性关系。此时很多人会从算法层面切入。

但是还有一种方式就是从硬件层面切入，使用合适的电路，让两个值之间是线性关系。