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 #endif   int    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输出波形频率。

如果波形一个周期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_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，但是实际用的时候无法满足单片机性能需求。单片机里面的算法，最好只有加减乘，不能有除。

离散傅里叶变换：

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;// really
	float 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);
		else
		printf("%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输出的正弦波的频率。

最后选择使用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);

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

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

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