基于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])实现频谱分析的功能，附带测试信号的生成代码。
- //基于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 {adc0,0,adc1,0,adc2,0,........}
- // 计算的频域数据也存储在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 NPT,const float df,float 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 ka,km;
- float alpha_m;
- float dkm;
-
- //加直流偏置
- for(i=0;i<NPT;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<NPT;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_len2048,x_X,0,1);
-
- //计算的频率、幅度、相位
- for (idx_harmonic=0;idx_harmonic<5;idx_harmonic++)
- {
- float magnitude_XH5_km,angle_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_k,Xreal_k)/PI*180;
- }
-
- //计算中心频率、幅度、相位
- //找到幅度极大谱线ka和临建次大谱线kb,取km=min(ka,kb)
- 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]<mag_neighborhood[ka+1-k0+harf_neighborhood])
- {
- km=ka;
- }
- else
- {
- km=ka-1;
- }
- alpha_m=mag_neighborhood[km-k0+harf_neighborhood]/mag_neighborhood[km+1-k0+harf_neighborhood];
- dkm=(4-3*alpha_m)/(1+alpha_m);
- magnitude_XH5_km=mag_neighborhood[km-k0+harf_neighborhood];
- angle_XH5_km=ang_neighborhood[km-k0+harf_neighborhood];
- }
- else
- {
- float Xreal_km,Ximag_km;
- float temp=f_harmonic[0]/df*(idx_harmonic+1);
- km=floor(temp);
- dkm=temp-km;
- 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.;
- 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.;
- Xreal_km/=256;//与(adc_value-2048)) << 8 对应
- Ximag_km/=256;//与(adc_value-2048)) << 8 对应
- magnitude_XH5_km=sqrt(Xreal_km*Xreal_km+Ximag_km*Ximag_km);//先不乘以2。实际上单边谱cn是双边谱Fn的2倍
- angle_XH5_km=atan2(Ximag_km,Xreal_km)/PI*180;
- }
- f_harmonic[idx_harmonic]=(km+dkm)*df;//中心频率
- //校正谐波参数
【必读】版权免责声明
1、本主题所有言论和内容纯属会员个人意见，与本论坛立场无关。2、本站对所发内容真实性、客观性、可用性不做任何保证也不负任何责任，网友之间仅出于学习目的进行交流。3、对提供的数字内容不拥有任何权利，其版权归原著者拥有。请勿将该数字内容进行商业交易、转载等行为，该内容只为学习所提供，使用后发生的一切问题与本站无关。 4、本网站不保证本站提供的下载资源的准确性、安全性和完整性;同时本网站也不承担用户因使用这些下载资源对自己和他人造成任何形式的损失或伤害。 5、本网站所有软件和资料均为网友推荐收集整理而来，仅供学习用途使用，请务必下载后两小时内删除，禁止商用。6、如有侵犯你版权的，请及时联系我们（电子邮箱1370723259@qq.com）指出，本站将立即改正。
|
|
