基于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）指出，本站將立即改正。