基于MATLAB的电力谐波分析，在有频谱泄露、栅栏效应

基于MATLAB FFT的电力谐波分析，输入时域信号采样值，在有频谱泄露、效应的前提下计算基波和谐波的频率、幅度、相位。根据论文《基于加汉宁窗的FFT高精度谐波检测改进算法》实现。函数[f_HA,A_HA,phi_HA]=harmonic_analysis(x,fs,N_HA)实现频谱分析的功能，主函数test_FFT产生仿真采样信号并调用频谱分析函数，计算电力质量指标。代码中有详细的注释。

1. function test_FFT
   
2. 
   
3. %信号采样的参数
   
4. N_HA=9;%待分析的谐波个数
   
5. fs=1024;%采样频率，大于谐波最高频率乘以2，至少N_HA*50*2或者N_HA*60*2
   
6. Nsample=2048;%1024,2048,4096
   
7. Nbit=12;%ADC采样的位数,取值8,10,12
   
8. kADC=280*sqrt(2)/2^(Nbit-1);%ADC的比例系数，允许最大电压为正负280*sqrt(2)
   
9. 
   
10. %电力质量数据缓存
    
11. Sizebuffer=10000;
    
12. U1buffer=zeros(1,Sizebuffer);%存储U1的缓存
    
13. f1buffer=zeros(1,Sizebuffer);%存储基波频率的缓存
    
14. harmonic_all_buffer=zeros(1,Sizebuffer);%总谐波比例
    
15. harmonic_even_buffer=zeros(1,Sizebuffer);%偶次谐波比例
    
16. harmonic_odd_buffer=zeros(1,Sizebuffer);%奇次谐波比例
    
17. 
    
18. for i=1:Sizebuffer
    
19.    
    
20.     %用采集卡采集电压数据
    
21.     %采样结果放到数组x中，x的数据点数为Nsample
    
22.     %TODO: 添加从采集卡中读取数据程序
    
23.    
    
24.    
    
25.    
    
26.     %软件测试时用matlab产生测试数据
    
27.     %电压信号的参数
    
28.     f1=50+randn()*0.3;%基波的频率，50,49.6,50.1
    
29.     %谐波（包含基波和高次谐波）的频率、幅度有效值（V）、初始相位（角度）
    
30.     f_HA_groundtruth=f1*(1:9);
    
31.     U1=220+randn()*5;
    
32.     Uharmonic=[5.3,3.2,2,9,1.1,8,1,0.5];%谐波的幅度
    
33.     Uharmonic(1:2:end)=Uharmonic(1:2:end)*rand()*1;%偶次谐波乘以系数，*0.5, *1, *3
    
34.     Uharmonic(2:2:end)=Uharmonic(2:2:end)*rand()*1;%奇次谐波乘以系数，*0.5, *1, *3
    
35.     A_HA_groundtruth=[U1,Uharmonic];
    
36.     phi_HA_groundtruth=[0,10,9,8,30,8,60,8,8];%有谐波信号
    
37.     %A_HA_groundtruth=[220,0,0,0,0,0,0,0,0];phi_HA_groundtruth=[0,0,0,0,0,0,0,0,0];%标准信号
    
38.     %生成采样数据
    
39.     t=(0:Nsample-1)/fs;%采样时刻数组
    
40.     x_contious=zeros(1,Nsample);
    
41.     for i_harmonic=1:numel(f_HA_groundtruth)
    
42.         x_contious=x_contious+A_HA_groundtruth(i_harmonic)*sqrt(2)*cos(f_HA_groundtruth(i_harmonic)*2*pi*t+phi_HA_groundtruth(i_harmonic)/180*pi);
    
43.     end
    
44.     %建模ADC采样的幅度离散化
    
45.     ADCvalue=round(x_contious/kADC);
    
46.     x=ADCvalue*kADC;
    
47.    
    
48.    
    
49.    
    
50. 
    
51.     %谐波分析
    
52.     [f_HA,A_HA,phi_HA]=harmonic_analysis(x,fs,N_HA);
    
53.     %比较计算结果
    
54.     %[f_HA;A_HA;phi_HA]
    
55.     %[f_HA_groundtruth;A_HA_groundtruth;phi_HA_groundtruth]
    
56.     U1buffer(i)=A_HA(1);
    
57.     f1buffer(i)=f_HA(1);
    
58.     harmonic_all_buffer(i)=norm(A_HA(2:end))/A_HA(1);
    
59.     harmonic_even_buffer(i)=norm(A_HA(2:2:end))/A_HA(1);
    
60.     harmonic_odd_buffer(i)=norm(A_HA(3:2:end))/A_HA(1);
    
61.    
    
62.     %结果分析
    
63.     figure(1),
    
64.     subplot(2,3,1)
    
65.     plot(U1buffer(1:i));
    
66.     ylabel('基波电压有效值(V)');
    
67.     subplot(2,3,2)
    
68.     plot(f1buffer(1:i));
    
69.     ylabel('基波频率(Hz)');
    
70.     subplot(2,3,4)
    
71.     plot(harmonic_all_buffer(1:i)*100);
    
72.     ylabel('总谐波比例(%)');
    
73.     subplot(2,3,5)
    
74.     plot(harmonic_even_buffer(1:i)*100);
    
75.     ylabel('偶次谐波比例(%)');
    
76.     subplot(2,3,6)
    
77.     plot(harmonic_odd_buffer(1:i)*100);
    
78.     ylabel('奇次谐波比例(%)');
    
79.    
    
80.     pause(0.5);%延时，使用数据采集卡时延时可以长些，模拟产生测试数据时减小延时
    
81. end
    
82. 
    
83. 
    
84. end
    
85. 
    
86. 
    
87. %电力谐波分析函数[f_HA,A_HA,phi_HA]=harmonic_analysis(x,fs,N_HA)
    
88. %根据论文基于加汉宁窗的FFT高精度谐波检测改进算法实现。
    
89. %在中心谱线估计中，针对两种极端情况进行近似计算，提高数值稳定性，避免得到错误结果。
    
90. %输入参数
    
91. %    x：数组长度为2^N，离散化后的电压采样值，单位为V。
    
92. %    fs:基波和各次谐波的频率。
    
93. %    N_HA:待分析的谐波个数。注意，f1*N_HA*2<fs。
    
94. %输出参数
    
95. %    f_HA：估计的谐波的频率，长度为N_HA。
    
96. %    A_HA：估计的谐波的幅度有效值（V)，长度为N_HA。
    
97. %    phi_HA：估计的谐波的角度，长度为N_HA。
    
98. function [f_HA,A_HA,phi_HA]=harmonic_analysis(x,fs,N_HA)
    
99. 
    
100. %检查采样频率大于信号最高频率的2倍
     
101. if N_HA*50*2>fs
     
102.     error('N_HA*50*2>fs');
     
103. end
     
104. 
     
105. f_HA=zeros(1,N_HA);
     
106. A_HA=zeros(1,N_HA);
     
107. 
     
108. %时域采样值加汉宁窗
     
109. Nsample=numel(x);
     
110. x=x.*hanning(Nsample).';
     
111. plot(x);
     
112. 
     
113. %计算加窗采样值的FFT
     
114. XH=fft(x);
     
115. plot(abs(XH))
     
116. plot(angle(XH)/pi*180)
     
117. 
     
118. %计算基波的频率、幅度、相位
     
119. %搜索主值谱线
     
120. df=fs/numel(x);%谱线间隔
     
121. k0=round(50/df)+1;%50Hz处的谱线位置
     
122. search_radius=round(15/df);%搜索半径，允许最大频偏15Hz
     
123. XH5_search=XH(k0-search_radius:k0+search_radius)/60 ...
     
124. -(XH(k0-1-search_radius:k0-1+search_radius)+XH(k0+1-search_radius:k0+1+search_radius))/90 ...
     
125. +(XH(k0-2-search_radius:k0-2+search_radius)+XH(k0+2-search_radius:k0+2+search_radius))/360;
     
126. mag_XH5_search=abs(XH5_search);
     
127. %ka,kb,km都是在mag_XH5_search中的下标
     
128. [~,ka]=max(mag_XH5_search);
     
129. if mag_XH5_search(ka-1)<mag_XH5_search(ka+1)
     
130. kb=ka+1;
     
131. else
     
132. kb=ka-1;
     
133. end
     
134. km=min(ka,kb);
     
135. %计算幅值比alpha_m
     
136. alpha_m=mag_XH5_search(km)/mag_XH5_search(km+1);
     
137. %计算偏移值dkm
     
138. dkm=(4-3*alpha_m)/(1+alpha_m);
     
139. %校正谐波参数
     
140. f_HA(1)=(km+k0-search_radius-1+dkm-1)*df;
     
141. if abs(dkm)<1e-3
     
142. A_HA(1)=abs((dkm^2-1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/pi*mag_XH5_search(km)/sqrt(2);%有效值
     
143. elseif abs(dkm-1)<1e-3
     
144. A_HA(1)=abs(dkm*(dkm+1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/pi*mag_XH5_search(km)/sqrt(2);%有效值
     
145. else
     
146. A_HA(1)=abs(dkm*(dkm^2-1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/sin(pi*dkm)*mag_XH5_search(km)/sqrt(2);%有效值
     
147. end
     
148. phi_HA(1)=angle(XH5_search(km))/pi*180-dkm*180;
     
149. 
     
150. %计算高次谐波的频率、幅度、相位
     
151. for order_HA=2:N_HA %谐波的序号
     
152. km=floor(f_HA(1)/df*order_HA+1);
     
153. dkm=f_HA(1)/df*order_HA+1-km;
     
154. 
     
155. XH5_km=XH(km)/60 ...
     
156. -(XH(km-1)+XH(km+1))/90 ...
     
157. +(XH(km-2)+XH(km+2))/360;
     
158. mag_XH5_km=abs(XH5_km);
     
159. f_HA(order_HA)=f_HA(1)*order_HA;%高次谐波的频率
     
160. if abs(dkm)<1e-3
     
161. A_HA(order_HA)=abs((dkm^2-1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/pi*mag_XH5_km/sqrt(2);%有效值
     
162. elseif abs(dkm-1)<1e-3
     
163. A_HA(order_HA)=abs(dkm*(dkm+1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/pi*mag_XH5_km/sqrt(2);%有效值
     
164. else
     
165. A_HA(order_HA)=abs(dkm*(dkm^2-1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/sin(pi*dkm)*mag_XH5_km/sqrt(2);%有效值
     
166. end
     
167. phi_HA(order_HA)=angle(XH5_km)/pi*180-dkm*180;
     
168. end
     
169. 
     
170. %限定相位角在[-180,180]
     
171. phi_HA(phi_HA<-180)=phi_HA(phi_HA<-180)+360;
     
172. 
     
173. end
     

复制代码

【必读】版权免责声明

1、本主题所有言论和内容纯属会员个人意见，与本论坛立场无关。2、本站对所发内容真实性、客观性、可用性不做任何保证也不负任何责任，网友之间仅出于学习目的进行交流。3、对提供的数字内容不拥有任何权利，其版权归原著者拥有。请勿将该数字内容进行商业交易、转载等行为，该内容只为学习所提供，使用后发生的一切问题与本站无关。 4、本网站不保证本站提供的下载资源的准确性、安全性和完整性;同时本网站也不承担用户因使用这些下载资源对自己和他人造成任何形式的损失或伤害。 5、本网站所有软件和资料均为网友推荐收集整理而来，仅供学习用途使用，请务必下载后两小时内删除，禁止商用。6、如有侵犯你版权的，请及时联系我们（电子邮箱1370723259@qq.com）指出，本站将立即改正。