基于STM32 DSP库的电力谐波分析根据ADC采样值计算基波谐波...

[复制链接]
查看2 | 回复0 | 2021-5-9 01:30:11 | 显示全部楼层 |阅读模式
基于STM32 DSP库函数的电力谐波分析,输入时域信号采样值,进行Q31 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])实现频谱分析的功能,附带测试信号的生成代码。
  1. //基于cfft的频谱分析,RAM空间复杂度8*NPT
  2. //频谱分析方法文献:基于加汉宁窗的FFT高精度谐波检测改进算法
  3. //reference: STM32的DSP库FFT函数 https://wenku.baidu.com/view/9cbc1a94eff9aef8951e061b.html
  4. //输入参数:
  5. //x_X: of size 2*NPT,  
  6. //        作为输入的ADC采样数据时,in the format of {adc0,0,adc1,0,adc2,0,........}
  7. //        计算的频域数据也存储在x_X中,  格式是{X[0]的实部,X[0]的实部,X[1]的实部,X[1]的虚部,X[2]的实部,X[2]的虚部,......}
  8. //df: 频率分辨率,等于采样频率除以采样点数
  9. //Magnitude_harmonic: 计算的基波及谐波的电压有效值
  10. //Angle_harmonic: 计算的基波及谐波的相位角
  11. //f_harmonic: 计算的基波及谐波的频率
  12. 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])
  13. {
  14.         int i;//时域信号的下标
  15.         int k;//频域信号的下标
  16.        
  17.         //谱线邻域分析用的变量
  18.         int k0;//谐波谱线的理论下标
  19.         const int harf_neighborhood=8;//正负15Hz邻域, 2/df
  20.         float mag_neighborhood[harf_neighborhood*2+1];
  21.         float ang_neighborhood[harf_neighborhood*2+1];
  22.         int idx_harmonic;//谐波的下标
  23.         float max_mag;
  24.         int ka,km;
  25.         float alpha_m;
  26.         float dkm;
  27.        
  28.         //加直流偏置
  29.         for(i=0;i<NPT;i++)
  30.         {
  31.                 //这里因为单片机的ADC只能测正的电压 所以需要前级加直流偏执
  32.                 //加入直流偏执后,需要在软件上减去2048即一半,达到负半周期测量的目的(需要根据具体情况来进行配置)
  33.                 //x[2*i] = ((signed short)(adc_buf[i]-2048)) << 16;
  34.                 //左移位数需要小于32-12,以避免定点FFT计算时溢出。
  35.                 x_X[2*i] = ((signed short)(x_X[2*i]-2048)) << 8;
  36.                 x_X[2*i+1]=0;
  37.         }
  38.        
  39.        
  40.         //数据加窗,汉宁窗
  41.         for(i=0;i<NPT;i++)
  42.                 x_X[2*i] = x_X[2*i]* 0.5*(1-cos(2*PI*i/(NPT+1)));
  43.        
  44.        
  45.         arm_cfft_q31(&arm_cfft_sR_q31_len2048,x_X,0,1);

  46.        
  47.         //计算的频率、幅度、相位
  48.         for (idx_harmonic=0;idx_harmonic<5;idx_harmonic++)
  49.         {
  50.                 float magnitude_XH5_km,angle_XH5_km;
  51.                 if (idx_harmonic==0)
  52.                 {
  53.                         //基波
  54.                         k0=(idx_harmonic+1)*50/df;//谐波的序号的理论值
  55.                         //计算谐波邻域内的幅度谱、相位谱
  56.                         for (k=k0-harf_neighborhood;k<=k0+harf_neighborhood;k++)
  57.                         {
  58.                                 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.;
  59.                                 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.;
  60.                                 Xreal_k/=256;//与(adc_value-2048)) << 8 对应
  61.                                 Ximag_k/=256;//与(adc_value-2048)) << 8 对应
  62.                                 mag_neighborhood[k-k0+harf_neighborhood]=sqrt(Xreal_k*Xreal_k+Ximag_k*Ximag_k);//先不乘以2。实际上单边谱cn是双边谱Fn的2倍
  63.                                 ang_neighborhood[k-k0+harf_neighborhood]=atan2(Ximag_k,Xreal_k)/PI*180;
  64.                         }
  65.                        
  66.                         //计算中心频率、幅度、相位
  67.                         //找到幅度极大谱线ka和临建次大谱线kb,取km=min(ka,kb)
  68.                         max_mag=0;
  69.                         for (k=k0-harf_neighborhood;k<=k0+harf_neighborhood;k++)
  70.                         {
  71.                                 if (mag_neighborhood[k-k0+harf_neighborhood]>max_mag)
  72.                                 {
  73.                                         ka=k;
  74.                                         max_mag=mag_neighborhood[k-k0+harf_neighborhood];
  75.                                 }
  76.                         }
  77.                         if (mag_neighborhood[ka-1-k0+harf_neighborhood]<mag_neighborhood[ka+1-k0+harf_neighborhood])
  78.                         {
  79.                                 km=ka;
  80.                         }
  81.                         else
  82.                         {
  83.                                 km=ka-1;
  84.                         }
  85.                         alpha_m=mag_neighborhood[km-k0+harf_neighborhood]/mag_neighborhood[km+1-k0+harf_neighborhood];
  86.                         dkm=(4-3*alpha_m)/(1+alpha_m);
  87.                         magnitude_XH5_km=mag_neighborhood[km-k0+harf_neighborhood];
  88.                         angle_XH5_km=ang_neighborhood[km-k0+harf_neighborhood];
  89.                 }
  90.                 else
  91.                 {
  92.                         float Xreal_km,Ximag_km;
  93.                         float temp=f_harmonic[0]/df*(idx_harmonic+1);
  94.                         km=floor(temp);
  95.                         dkm=temp-km;
  96.                         Xreal_km=x_X[2*km]/60.-(x_X[2*(km+1)]+x_X[2*(km-1)])/90.+(x_X[2*(km+2)]+x_X[2*(km-2)])/360.;
  97.                         Ximag_km=x_X[2*km+1]/60.-(x_X[2*(km+1)+1]+x_X[2*(km-1)+1])/90.+(x_X[2*(km+2)+1]+x_X[2*(km-2)+1])/360.;
  98.                         Xreal_km/=256;//与(adc_value-2048)) << 8 对应
  99.                         Ximag_km/=256;//与(adc_value-2048)) << 8 对应
  100.                         magnitude_XH5_km=sqrt(Xreal_km*Xreal_km+Ximag_km*Ximag_km);//先不乘以2。实际上单边谱cn是双边谱Fn的2倍
  101.                         angle_XH5_km=atan2(Ximag_km,Xreal_km)/PI*180;
  102.                 }
  103.                 f_harmonic[idx_harmonic]=(km+dkm)*df;//中心频率
  104.                 //校正谐波参数
复制代码

剩余代码:在下方已被隐藏,消耗积分可见
****查看此区域内容****    需支付 2 工控币后可查看立即购买
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则