基于MATLAB的電力諧波分析，在有頻譜泄露、柵欄效應 - ELEOK

function test_FFT
%信號采樣的參數
N_HA=9;%待分析的諧波個數
fs=1024;%采樣頻率，大于諧波最高頻率乘以2，至少N_HA*50*2或者N_HA*60*2
Nsample=2048;%1024,2048,4096
Nbit=12;%ADC采樣的位數,取值8,10,12
kADC=280*sqrt(2)/2^(Nbit-1);%ADC的比例系數，允許最大電壓為正負280*sqrt(2)
%電力質量數據緩存
Sizebuffer=10000;
U1buffer=zeros(1,Sizebuffer);%存儲U1的緩存
f1buffer=zeros(1,Sizebuffer);%存儲基波頻率的緩存
harmonic_all_buffer=zeros(1,Sizebuffer);%總諧波比例
harmonic_even_buffer=zeros(1,Sizebuffer);%偶次諧波比例
harmonic_odd_buffer=zeros(1,Sizebuffer);%奇次諧波比例
for i=1:Sizebuffer
%用采集卡采集電壓數據
%采樣結果放到數組x中，x的數據點數為Nsample
%TODO: 添加從采集卡中讀取數據程序
%軟件測試時用matlab產生測試數據
%電壓信號的參數
f1=50+randn()*0.3;%基波的頻率，50,49.6,50.1
%諧波（包含基波和高次諧波）的頻率、幅度有效值（V）、初始相位（角度）
f_HA_groundtruth=f1*(1:9);
U1=220+randn()*5;
Uharmonic=[5.3,3.2,2,9,1.1,8,1,0.5];%諧波的幅度
Uharmonic(1:2:end)=Uharmonic(1:2:end)*rand()*1;%偶次諧波乘以系數，*0.5, *1, *3
Uharmonic(2:2:end)=Uharmonic(2:2:end)*rand()*1;%奇次諧波乘以系數，*0.5, *1, *3
A_HA_groundtruth=[U1,Uharmonic];
phi_HA_groundtruth=[0,10,9,8,30,8,60,8,8];%有諧波信號
%A_HA_groundtruth=[220,0,0,0,0,0,0,0,0];phi_HA_groundtruth=[0,0,0,0,0,0,0,0,0];%標準信號
%生成采樣數據
t=(0:Nsample-1)/fs;%采樣時刻數組
x_contious=zeros(1,Nsample);
for i_harmonic=1:numel(f_HA_groundtruth)
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);
end
%建模ADC采樣的幅度離散化
ADCvalue=round(x_contious/kADC);
x=ADCvalue*kADC;
%諧波分析
[f_HA,A_HA,phi_HA]=harmonic_analysis(x,fs,N_HA);
%比較計算結果
%[f_HA;A_HA;phi_HA]
%[f_HA_groundtruth;A_HA_groundtruth;phi_HA_groundtruth]
U1buffer(i)=A_HA(1);
f1buffer(i)=f_HA(1);
harmonic_all_buffer(i)=norm(A_HA(2:end))/A_HA(1);
harmonic_even_buffer(i)=norm(A_HA(2:2:end))/A_HA(1);
harmonic_odd_buffer(i)=norm(A_HA(3:2:end))/A_HA(1);
%結果分析
figure(1),
subplot(2,3,1)
plot(U1buffer(1:i));
ylabel('基波電壓有效值(V)');
subplot(2,3,2)
plot(f1buffer(1:i));
ylabel('基波頻率(Hz)');
subplot(2,3,4)
plot(harmonic_all_buffer(1:i)*100);
ylabel('總諧波比例(%)');
subplot(2,3,5)
plot(harmonic_even_buffer(1:i)*100);
ylabel('偶次諧波比例(%)');
subplot(2,3,6)
plot(harmonic_odd_buffer(1:i)*100);
ylabel('奇次諧波比例(%)');
pause(0.5);%延時，使用數據采集卡時延時可以長些，模擬產生測試數據時減小延時
end
end
%電力諧波分析函數[f_HA,A_HA,phi_HA]=harmonic_analysis(x,fs,N_HA)
%根據論文基于加漢寧窗的FFT高精度諧波檢測改進算法實現。
%在中心譜線估計中，針對兩種極端情況進行近似計算，提高數值穩定性，避免得到錯誤結果。
%輸入參數
% x：數組長度為2^N，離散化后的電壓采樣值，單位為V。
% fs:基波和各次諧波的頻率。
% N_HA:待分析的諧波個數。注意，f1*N_HA*2<fs。
%輸出參數
% f_HA：估計的諧波的頻率，長度為N_HA。
% A_HA：估計的諧波的幅度有效值（V)，長度為N_HA。
% phi_HA：估計的諧波的角度，長度為N_HA。
function [f_HA,A_HA,phi_HA]=harmonic_analysis(x,fs,N_HA)
%檢查采樣頻率大于信號最高頻率的2倍
if N_HA*50*2>fs
error('N_HA*50*2>fs');
end
f_HA=zeros(1,N_HA);
A_HA=zeros(1,N_HA);
%時域采樣值加漢寧窗
Nsample=numel(x);
x=x.*hanning(Nsample).';
plot(x);
%計算加窗采樣值的FFT
XH=fft(x);
plot(abs(XH))
plot(angle(XH)/pi*180)
%計算基波的頻率、幅度、相位
%搜索主值譜線
df=fs/numel(x);%譜線間隔
k0=round(50/df)+1;%50Hz處的譜線位置
search_radius=round(15/df);%搜索半徑，允許最大頻偏15Hz
XH5_search=XH(k0-search_radius:k0+search_radius)/60 ...
-(XH(k0-1-search_radius:k0-1+search_radius)+XH(k0+1-search_radius:k0+1+search_radius))/90 ...
+(XH(k0-2-search_radius:k0-2+search_radius)+XH(k0+2-search_radius:k0+2+search_radius))/360;
mag_XH5_search=abs(XH5_search);
%ka,kb,km都是在mag_XH5_search中的下標
[~,ka]=max(mag_XH5_search);
if mag_XH5_search(ka-1)<mag_XH5_search(ka+1)
kb=ka+1;
else
kb=ka-1;
end
km=min(ka,kb);
%計算幅值比alpha_m
alpha_m=mag_XH5_search(km)/mag_XH5_search(km+1);
%計算偏移值dkm
dkm=(4-3*alpha_m)/(1+alpha_m);
%校正諧波參數
f_HA(1)=(km+k0-search_radius-1+dkm-1)*df;
if abs(dkm)<1e-3
A_HA(1)=abs((dkm^2-1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/pi*mag_XH5_search(km)/sqrt(2);%有效值
elseif abs(dkm-1)<1e-3
A_HA(1)=abs(dkm*(dkm+1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/pi*mag_XH5_search(km)/sqrt(2);%有效值
else
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);%有效值
end
phi_HA(1)=angle(XH5_search(km))/pi*180-dkm*180;
%計算高次諧波的頻率、幅度、相位
for order_HA=2:N_HA %諧波的序號
km=floor(f_HA(1)/df*order_HA+1);
dkm=f_HA(1)/df*order_HA+1-km;
XH5_km=XH(km)/60 ...
-(XH(km-1)+XH(km+1))/90 ...
+(XH(km-2)+XH(km+2))/360;
mag_XH5_km=abs(XH5_km);
f_HA(order_HA)=f_HA(1)*order_HA;%高次諧波的頻率
if abs(dkm)<1e-3
A_HA(order_HA)=abs((dkm^2-1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/pi*mag_XH5_km/sqrt(2);%有效值
elseif abs(dkm-1)<1e-3
A_HA(order_HA)=abs(dkm*(dkm+1)*(dkm^2-4)*(dkm^2-9))*4*pi/Nsample/pi*mag_XH5_km/sqrt(2);%有效值
else
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);%有效值
end
phi_HA(order_HA)=angle(XH5_km)/pi*180-dkm*180;
end
%限定相位角在[-180,180]
phi_HA(phi_HA<-180)=phi_HA(phi_HA<-180)+360;
end

復制代碼

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

ELEOK