基于STM32 DSP库的电力谐波分析，根据ADC采样值计算基波谐波的频率幅度相位

基于STM32 DSP库函数的电力谐波分析，输入时域信号采样值，进行Q31 FFT计算，在有频谱泄露、栅栏效应的前提下计算基波和谐波的频率、幅度、相位。根据论文《基于加汉宁窗的FFT高精度谐波检测改进算法》实现。函数void spectrum_with_cfft(q31_t * x_X, const int NPT,const float df,float Magnitude_harmonic[5],float Angle_harmonic[5],float f_harmonic[5])实现频谱分析的功能，附带测试信号的生成代码

//基于cfft的频谱分析，RAM空间复杂度8*NPT
//频谱分析方法文献：基于加汉宁窗的FFT高精度谐波检测改进算法
//reference: STM32的DSP库FFT函数 https://wenku.baidu.com/view/9cbc1a94eff9aef8951e061b.html
//输入参数：
//x_X: of size 2*NPT  
//        作为输入的ADC采样数据时，in the format of {adc00adc10adc20........}
//        计算的频域数据也存储在x_X中  格式是{X[0]的实部，X[0]的实部，X[1]的实部，X[1]的虚部，X[2]的实部，X[2]的虚部，......}
//df: 频率分辨率，等于采样频率除以采样点数
//Magnitude_harmonic: 计算的基波及谐波的电压有效值
//Angle_harmonic: 计算的基波及谐波的相位角
//f_harmonic: 计算的基波及谐波的频率
void spectrum_with_cfft（q31_t * x_X const int NPTconst float dffloat Magnitude_harmonic[5]float Angle_harmonic[5]float f_harmonic[5]）
{
	int i;//时域信号的下标
	int k;//频域信号的下标
	
	//谱线邻域分析用的变量
	int k0;//谐波谱线的理论下标
	const int harf_neighborhood=8;//正负15Hz邻域 2/df
	float mag_neighborhood[harf_neighborhood*2+1];
	float ang_neighborhood[harf_neighborhood*2+1];
	int idx_harmonic;//谐波的下标
	float max_mag;
	int kakm;
	float alpha_m;
	float dkm;
	
	//加直流偏置
	for（i=0;i	{
		//这里因为单片机的ADC只能测正的电压 所以需要前级加直流偏执
		//加入直流偏执后，需要在软件上减去2048即一半，达到负半周期测量的目的（需要根据具体情况来进行配置）
		//x[2*i] = （（signed short）（adc_buf[i]-2048）） << 16;
		//左移位数需要小于32-12，以避免定点FFT计算时溢出。
		x_X[2*i] = （（signed short）（x_X[2*i]-2048）） << 8;
		x_X[2*i+1]=0;
	}
	
	
	//数据加窗，汉宁窗
	for（i=0;i		x_X[2*i] = x_X[2*i]* 0.5*（1-cos（2*PI*i/（NPT+1）））;
	
	
	arm_cfft_q31（&arm_cfft_sR_q31_len2048x_X01）;

	
	//计算的频率、幅度、相位
	for （idx_harmonic=0;idx_harmonic<5;idx_harmonic++）
	{
		float magnitude_XH5_kmangle_XH5_km;
		if （idx_harmonic==0）
		{
			//基波
			k0=（idx_harmonic+1）*50/df;//谐波的序号的理论值
			//计算谐波邻域内的幅度谱、相位谱
			for （k=k0-harf_neighborhood;k<=k0+harf_neighborhood;k++）
			{
				float Xreal_k=x_X[2*k]/60.-（x_X[2*（k+1）]+x_X[2*（k-1）]）/90.+（x_X[2*（k+2）]+x_X[2*（k-2）]）/360.;
				float Ximag_k=x_X[2*k+1]/60.-（x_X[2*（k+1）+1]+x_X[2*（k-1）+1]）/90.+（x_X[2*（k+2）+1]+x_X[2*（k-2）+1]）/360.;
				Xreal_k/=256;//与（adc_value-2048）） << 8 对应
				Ximag_k/=256;//与（adc_value-2048）） << 8 对应
				mag_neighborhood[k-k0+harf_neighborhood]=sqrt（Xreal_k*Xreal_k+Ximag_k*Ximag_k）;//先不乘以2。实际上单边谱cn是双边谱Fn的2倍
				ang_neighborhood[k-k0+harf_neighborhood]=atan2（Ximag_kXreal_k）/PI*180;
			}
			
			//计算中心频率、幅度、相位
			//找到幅度极大谱线ka和临建次大谱线kb取km=min（kakb）
			max_mag=0;
			for （k=k0-harf_neighborhood;k<=k0+harf_neighborhood;k++）
			{
				if （mag_neighborhood[k-k0+harf_neighborhood]>max_mag）
				{
					ka=k;
					max_mag=mag_neighborhood[k-k0+harf_neighborhood];
				}
			}
			if （mag_neighborhood[ka-1-k0+harf_neighborhood]			{
				km=ka;
			}
			else
			{
				km=ka-1;
			}
			alpha_m=mag_neighborhood[km-k0+harf_neighborhood]/mag_n