国产免费AV|泡泡玛特欧洲总部将设在伦敦|中文天堂网www新版资源在线|一本久道综合在线中文|国精产品一二三产区的使用方法|香蕉鱼在线观看|www.27eee

ELEOK

標題: 基于STM32 DSP庫的電力諧波分析根據ADC采樣值計算基波諧波... [打印本頁]
作者: bosa    時間: 2021-5-9 01:30
標題: 基于STM32 DSP庫的電力諧波分析根據ADC采樣值計算基波諧波...
基于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);

  47.         
  48.         //計算的頻率、幅度、相位
  49.         for (idx_harmonic=0;idx_harmonic<5;idx_harmonic++)
  50.         {
  51.                 float magnitude_XH5_km,angle_XH5_km;
  52.                 if (idx_harmonic==0)
  53.                 {
  54.                         //基波
  55.                         k0=(idx_harmonic+1)*50/df;//諧波的序號的理論值
  56.                         //計算諧波鄰域內的幅度譜、相位譜
  57.                         for (k=k0-harf_neighborhood;k<=k0+harf_neighborhood;k++)
  58.                         {
  59.                                 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.;
  60.                                 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.;
  61.                                 Xreal_k/=256;//與(adc_value-2048)) << 8 對應
  62.                                 Ximag_k/=256;//與(adc_value-2048)) << 8 對應
  63.                                 mag_neighborhood[k-k0+harf_neighborhood]=sqrt(Xreal_k*Xreal_k+Ximag_k*Ximag_k);//先不乘以2。實際上單邊譜cn是雙邊譜Fn的2倍
  64.                                 ang_neighborhood[k-k0+harf_neighborhood]=atan2(Ximag_k,Xreal_k)/PI*180;
  65.                         }
  66.                         
  67.                         //計算中心頻率、幅度、相位
  68.                         //找到幅度極大譜線ka和臨建次大譜線kb,取km=min(ka,kb)
  69.                         max_mag=0;
  70.                         for (k=k0-harf_neighborhood;k<=k0+harf_neighborhood;k++)
  71.                         {
  72.                                 if (mag_neighborhood[k-k0+harf_neighborhood]>max_mag)
  73.                                 {
  74.                                         ka=k;
  75.                                         max_mag=mag_neighborhood[k-k0+harf_neighborhood];
  76.                                 }
  77.                         }
  78.                         if (mag_neighborhood[ka-1-k0+harf_neighborhood]<mag_neighborhood[ka+1-k0+harf_neighborhood])
  79.                         {
  80.                                 km=ka;
  81.                         }
  82.                         else
  83.                         {
  84.                                 km=ka-1;
  85.                         }
  86.                         alpha_m=mag_neighborhood[km-k0+harf_neighborhood]/mag_neighborhood[km+1-k0+harf_neighborhood];
  87.                         dkm=(4-3*alpha_m)/(1+alpha_m);
  88.                         magnitude_XH5_km=mag_neighborhood[km-k0+harf_neighborhood];
  89.                         angle_XH5_km=ang_neighborhood[km-k0+harf_neighborhood];
  90.                 }
  91.                 else
  92.                 {
  93.                         float Xreal_km,Ximag_km;
  94.                         float temp=f_harmonic[0]/df*(idx_harmonic+1);
  95.                         km=floor(temp);
  96.                         dkm=temp-km;
  97.                         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.;
  98.                         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.;
  99.                         Xreal_km/=256;//與(adc_value-2048)) << 8 對應
  100.                         Ximag_km/=256;//與(adc_value-2048)) << 8 對應
  101.                         magnitude_XH5_km=sqrt(Xreal_km*Xreal_km+Ximag_km*Ximag_km);//先不乘以2。實際上單邊譜cn是雙邊譜Fn的2倍
  102.                         angle_XH5_km=atan2(Ximag_km,Xreal_km)/PI*180;
  103.                 }
  104.                 f_harmonic[idx_harmonic]=(km+dkm)*df;//中心頻率
  105.                 //校正諧波參數
復制代碼

作者: googzdz    時間: 2023-5-18 15:16
感謝樓主發布
作者: bison    時間: 2024-3-26 14:28
感謝分析



歡迎光臨 ELEOK (http://m.afoofa.cn/) Powered by Discuz! X5.0