Code covered by the BSD License

# IEC 61000-4-7 for Harmonics Calculations

by

### Aubai (view profile)

This file gives you the ability to calculated Harmonics

IEC_Harmonics_final(f,Vn,Inom,U1,U2,U3,I1,I2,I3,U1_time,U2_time,U3_time,I1_time,I2_time,I3_time,figures)
```function [THD_U,THD_I,K_U,K_I,MaxHarm_I,MaxIHarm_I,MaxHFHarm_I,MaxHarm_U,MaxIHarm_U,MaxHFHarm_U,Order_Harm,Order_IHarm,Order_HFHarm] = IEC_Harmonics_final(f,Vn,Inom,U1,U2,U3,I1,I2,I3,U1_time,U2_time,U3_time,I1_time,I2_time,I3_time,figures)
% This function will calculate the Harmonics spectrum values according to
% the IEC 61000-4-7 Norm which spesify the need for 10 min Voltage and
% current measured signals of the system under study and calulate the
% Harmonics mean values every 10 Periods of the measured signal (In this case 50 Hz signal)
% This Fucntion inputs are:
%           1- The nominal values of the measured signal: Frequency,
%           Voltage, and the current. For example f = 50%Hz
%                                                 Vn = 400%V (Phase to Phase Voltage of the grid)
%                                                 Inom = 475%A (The nominal
%                                                 current of the system under measurement)
%           2- The measured Voltages and current signals with there time
%           measurement signals (U1,U2,U3,I1,I2,I3,U1_time,U2_time,U3_time,I1_time,I2_time,I3_time)
%           3- figures is an input of 1 or 0 in order to activate showing the figures of the different  harmonics spectrums figures calculated in this function.
%PS: This function assume that you have the same sampling rate of the
%different measured signals and that all measured signals has the same
%length.
Fn = f;
Un = Vn;
In = Inom;
Delta_Time_Orginal = I1_time(2)-I1_time(1);
F_Sample = abs(1/Delta_Time_Orginal);
%------------------------------------------------------------------------------------------------------
%---------------------- The Harmonics Spectrum of the complet input signal ----------------------------
%------------------------------------------------------------------------------------------------------
Tp = length(I3_time);
T_S = abs(Tp*Delta_Time_Orginal);
[nyquist_number_I1,Bezugsfrequenz_I1,Signal_amp_I1] = Fast_Harmonics (I1,F_Sample,T_S);
[nyquist_number_I2,Bezugsfrequenz_I2,Signal_amp_I2] = Fast_Harmonics (I2,F_Sample,T_S);
[nyquist_number_I3,Bezugsfrequenz_I3,Signal_amp_I3] = Fast_Harmonics (I3,F_Sample,T_S);
[nyquist_number_U1,Bezugsfrequenz_U1,Signal_amp_U1] = Fast_Harmonics (U1,F_Sample,T_S);
[nyquist_number_U2,Bezugsfrequenz_U2,Signal_amp_U2] = Fast_Harmonics (U2,F_Sample,T_S);
[nyquist_number_U3,Bezugsfrequenz_U3,Signal_amp_U3] = Fast_Harmonics (U3,F_Sample,T_S);
if figures
figure;
subplot(2,1,1);
plot([0:nyquist_number_I1]*Bezugsfrequenz_I1,Signal_amp_I1);ylabel('amplitude Currents [A]');title('Harmonic spectrom');hold on
plot([0:nyquist_number_I2]*Bezugsfrequenz_I2,Signal_amp_I2,'-.r');hold on
plot([0:nyquist_number_I3]*Bezugsfrequenz_I3,Signal_amp_I3,'--k');hold off
subplot(2,1,2);
plot([0:nyquist_number_U1]*Bezugsfrequenz_U1,Signal_amp_U1);xlabel('frequency [Hz]'),ylabel('amplitude Voltages [V]');hold on
plot([0:nyquist_number_U2]*Bezugsfrequenz_U2,Signal_amp_U2,'-.r');hold on
plot([0:nyquist_number_U3]*Bezugsfrequenz_U3,Signal_amp_U3,'.k','MarkerSize',0.1);hold off
end
%**********************************************************************
%---------------------- Creating a grouping matrix --------------------
%**********************************************************************
Sampling_Time = Delta_Time_Orginal;
T_10_period = 10/Fn;
F_Sample_Orginal = 1/Sampling_Time;
if abs(F_Sample_Orginal) < 20000
F_Sample_Complete_Data = 20000;
[p, q] = rat(abs(Delta_Time_Orginal/(1/F_Sample_Complete_Data)));
I1_Resampled = resample(I1,p,q);
I2_Resampled = resample(I2,p,q);
I3_Resampled = resample(I3,p,q);
U1_Resampled = resample(U1,p,q);
U2_Resampled = resample(U2,p,q);
U3_Resampled = resample(U3,p,q);
t_P = I1_time(1):((q/p)/F_Sample_Orginal):(I1_time(end)+((q/p)/F_Sample_Orginal));
tt = length(t_P);
if length(t_P) == length(I1_Resampled)
I1_Resampled_time = t_P;
I2_Resampled_time = t_P;
I3_Resampled_time = t_P;
U1_Resampled_time  = t_P;
U2_Resampled_time  = t_P;
U3_Resampled_time  = t_P;
else if length(t_P) < length(I1_Resampled)
for kk = 1:abs(length(I1_Resampled)-length(t_P))
t_P(tt+kk) = ((q/p)/F_Sample_Orginal)+t_P(tt+kk-1);
end
I1_Resampled_time = t_P;
I2_Resampled_time = t_P;
I3_Resampled_time = t_P;
U1_Resampled_time  = t_P;
U2_Resampled_time  = t_P;
U3_Resampled_time  = t_P;
elseif length(t_P) > length(I1_Resampled)
t_P = t_P(1:length(I1_Resample));
I1_Resampled_time = t_P;
I2_Resampled_time = t_P;
I3_Resampled_time = t_P;
U1_Resampled_time  = t_P;
U2_Resampled_time  = t_P;
U3_Resampled_time  = t_P;
end
end
Delta_Time_Complete_Data = 1/F_Sample_Complete_Data;
Sampling_Time_New = Delta_Time_Complete_Data;
Sampling_frq_New = (1/Sampling_Time_New)/1000;
Sample_New = T_10_period*Sampling_frq_New*1000;
ii= 1;
Points = 2^ii;
while Points < Sample_New
Points = 2^ii;
ii = ii+1;
end
else if abs(F_Sample_Orginal) > 20000
F_Sample_Complete_Data = 20000;
[p, q] = rat(abs(Delta_Time_Orginal/(1/F_Sample_Complete_Data)));
I1_Resampled = resample(I1,p,q);
I2_Resampled = resample(I2,p,q);
I3_Resampled = resample(I3,p,q);
U1_Resampled = resample(U1,p,q);
U2_Resampled = resample(U2,p,q);
U3_Resampled = resample(U3,p,q);
t_Q = I1_time(1):((q/p)/F_Sample):I1_time(end);
I1_Resampled_time = t_Q;
I2_Resampled_time = t_Q;
I3_Resampled_time = t_Q;
U1_Resampled_time  = t_Q;
U2_Resampled_time  = t_Q;
U3_Resampled_time  = t_Q;
Delta_Time_Complete_Data = 1/F_Sample_Complete_Data;
%Sampling_Time_New = Delta_Time_Complete_Data;
Sampling_frq_New = (1/Sampling_Time_Complete_Data)/1000;
Sample_New = T_10_period*Sampling_frq_New*1000;
ii= 1;
Points = 2^ii;
while Points < Sample_New
Points = 2^ii;
ii = ii+1;
end
else
I1_Resampled = I1;
I2_Resampled = I2;
I3_Resampled = I3;
U1_Resampled = U1;
U2_Resampled = U2;
U3_Resampled = U3;
I1_Resampled_time = I1_time;
I2_Resampled_time = I2_time;
I3_Resampled_time = I3_time;
U1_Resampled_time  = U1_time;
U2_Resampled_time  = U2_time;
U3_Resampled_time  = U3_time;
F_Sample_Complete_Data = 20000;
Delta_Time_Complete_Data = I1_Resampled_time(2)-I1_Resampled_time(1);
Sampling_frq = 1/Delta_Time_Complete_Data;
T_10_period = 10/Fn;
Sample = T_10_period*Sampling_frq;
ii= 1;
Points = 2^ii;
while Points < Sample
Points = 2^ii;
ii = ii+1;
end
end
end
%----------------------------------------------------------------------
%--------------- Zeros Results Array Definition -----------------------
%----------------------------------------------------------------------
%--------------- Try --------------------
LOWPASS_CUTOFF = 100;% Filter cutoff frequency is 100 Hz
LOWPASS_ORDER = 4;% Filter order is 4
[b_bw, a_bw] = butter(LOWPASS_ORDER, LOWPASS_CUTOFF / (F_Sample_Complete_Data / 2), 'low');
Filter = filter(b_bw, a_bw, U1_Resampled);
%--------------- End --------------------
null_L1 = find(diff(Filter < 0));
null_L1_x = Zeros_finding(Filter,I1_Resampled_time');
Tp = length(U1_Resampled_time);
len = length(null_L1_x);
%len_1_period = (1/Fn)/Delta_Time_Complete_Data;
%len_10_period = (10*(1/Fn))/Delta_Time_Complete_Data;
Results_U_H1 = zeros(50,length(0:20:len-30));
Results_U_H2 = zeros(50,length(0:20:len-30));
Results_U_H3 = zeros(50,length(0:20:len-30));
Results_U_IH1 = zeros(39,length(0:20:len-30));
Results_U_IH2 = zeros(39,length(0:20:len-30));
Results_U_IH3 = zeros(39,length(0:20:len-30));
Results_U_HF1 = zeros(35,length(0:20:len-30));
Results_U_HF2 = zeros(35,length(0:20:len-30));
Results_U_HF3 = zeros(35,length(0:20:len-30));
Results_P = zeros(length(0:20:len-30));
Results_Q = zeros(length(0:20:len-30));
%----------------------------------------------------------------------
%--------------------------- End --------------------------------------
%----------------------------------------------------------------------
%**********************************************************************
%----------------------------------END---------------------------------
%**********************************************************************
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%**********************************************************************
%--------------------------- Gruoping ---------------------------------
%**********************************************************************
j = 0;
jj = 1;
while j < len-30 % t < t2+0.1
data_I1 = I1_Resampled(null_L1(j+3):null_L1(j+23)+1);
data_I2 = I2_Resampled(null_L1(j+3):null_L1(j+23)+1);
data_I3 = I3_Resampled(null_L1(j+3):null_L1(j+23)+1);
data_U1 = U1_Resampled(null_L1(j+3):null_L1(j+23)+1);
data_U2 = U2_Resampled(null_L1(j+3):null_L1(j+23)+1);
data_U3 = U3_Resampled(null_L1(j+3):null_L1(j+23)+1);
data_I1_time = I1_Resampled_time(null_L1(j+3):null_L1(j+23)+1);
%------------------- Cut End --------------------------------------
%------------------------------------------------------------------
%--------- Data Resampling with the new Delta_10_Periods ----------
%------------------------------------------------------------------
delta = (null_L1_x(j+23)-null_L1_x(j+3))/Points;
%Rampe = null_L1_x(j+3):delta:(((Points-1)*delta)+null_L1_x(j+3));
iii=2;
Rampe = zeros();
Rampe(1) = null_L1_x(j+3);
while iii <= Points
Rampe(iii) = delta+Rampe(iii-1);
iii = iii+1;
end
F_Sample_10_Periods = 1/delta;
if abs(F_Sample_10_Periods) < F_Sample_Complete_Data
F_Sample_New = F_Sample_10_Periods;
[p, q] = rat(abs(Delta_Time/(1/F_Sample_New)));
I1_Resampled_10_Period = resample(data_I1,p,q);
I2_Resampled_10_Period = resample(data_I2,p,q);
I3_Resampled_10_Period = resample(data_I3,p,q);
U1_Resampled_10_Period = resample(data_U1,p,q);
U2_Resampled_10_Period = resample(data_U2,p,q);
U3_Resampled_10_Period = resample(data_U3,p,q);
t_P = data_I1_time(1):(1/F_Sample_New):(data_I1_time(end)+(1/F_Sample_New));
I1_Resampled_time_10_Period = t_P;
I2_Resampled_time_10_Period = t_P;
I3_Resampled_time_10_Period = t_P;
U1_Resampled_time_10_Period  = t_P;
U2_Resampled_time_10_Period  = t_P;
U3_Resampled_time_10_Period  = t_P;
else if abs(F_Sample_10_Periods) > F_Sample_Complete_Data
F_Sample_New = F_Sample_10_Periods;
[p, q] = rat(abs(Delta_Time_Complete_Data/(1/F_Sample_New)));
I1_Resampled_10_Period = resample(data_I1,p,q);
I2_Resampled_10_Period = resample(data_I2,p,q);
I3_Resampled_10_Period = resample(data_I3,p,q);
U1_Resampled_10_Period = resample(data_U1,p,q);
U2_Resampled_10_Period = resample(data_U2,p,q);
U3_Resampled_10_Period = resample(data_U3,p,q);
for kk=1:length(I1_Resampled_10_Period)-1
I1_Resample_10_Period_Final(kk) =  ((I1_Resampled_10_Period(kk)+I1_Resampled_10_Period(kk+1))/2);
I2_Resample_10_Period_Final(kk) =  ((I2_Resampled_10_Period(kk)+I2_Resampled_10_Period(kk+1))/2);
I3_Resample_10_Period_Final(kk) =  ((I3_Resampled_10_Period(kk)+I3_Resampled_10_Period(kk+1))/2);
U1_Resample_10_Period_Final(kk) =  ((U1_Resampled_10_Period(kk)+U1_Resampled_10_Period(kk+1))/2);
U2_Resample_10_Period_Final(kk) =  ((U2_Resampled_10_Period(kk)+U2_Resampled_10_Period(kk+1))/2);
U3_Resample_10_Period_Final(kk) =  ((U3_Resampled_10_Period(kk)+U3_Resampled_10_Period(kk+1))/2);
end
if length(I1_Resample_10_Period_Final) > length(Rampe)
I1_Resample_10_Period_Final = I1_Resample_10_Period_Final(1:length(Rampe));
I2_Resample_10_Period_Final = I2_Resample_10_Period_Final(1:length(Rampe));
I3_Resample_10_Period_Final = I3_Resample_10_Period_Final(1:length(Rampe));
else if length(I1_Resample_10_Period_Final) < length(Rampe)
error('Please check the Resampled data something wrong there :(')
end
end
t_Q = Rampe(1):((q/p)/F_Sample_Complete_Data):Rampe(end);
if length(t_Q) == length(I1_Resampled_10_Period)
I1_Resampled_time_10_Period = t_Q;
I2_Resampled_time_10_Period = t_Q;
I3_Resampled_time_10_Period = t_Q;
U1_Resampled_time_10_Period  = t_Q;
U2_Resampled_time_10_Period  = t_Q;
U3_Resampled_time_10_Period  = t_Q;
else if length(t_Q) > length(I1_Resampled_10_Period)
I1_Resampled_time_10_Period = t_Q(1:length(I1_Resampled_10_Period));
I2_Resampled_time_10_Period = t_Q(1:length(I2_Resampled_10_Period));
I3_Resampled_time_10_Period = t_Q(1:length(I3_Resampled_10_Period));
U1_Resampled_time_10_Period = t_Q(1:length(U1_Resampled_10_Period));
U2_Resampled_time_10_Period = t_Q(1:length(U2_Resampled_10_Period));
U3_Resampled_time_10_Period = t_Q(1:length(U3_Resampled_10_Period));
else if length(t_Q) < length(I1_Resampled_10_Period)
I1_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
I2_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
I3_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
U1_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
U2_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
U3_Resampled_time_10_Period = [t_Q t_Q(end)+((q/p)/F_Sample_Complete_Data)];
end
end
end
else
error('Please check the Delta Value of the 10 Period signal')
end
end
%------------------------------------------------------------------
%---------------- Data Re-sampling End ----------------------------
%------------------------------------------------------------------
%----------------Spectrum Calculation------------------------------
%------------------------------------------------------------------
Tp_10_Period_Resampled = length(I1_Resampled_time_10_Period);
Delta_Time_10_Period_Resampled = I1_Resampled_time_10_Period(2)-I1_Resampled_time_10_Period(3);
F_Sample_10_Period_Resampled = abs(1/Delta_Time_10_Period_Resampled);
T_S_10_Period_Resampled = abs(Tp*Delta_Time_10_Period_Resampled);
%----------------------------------------------------------------------
[nyquist_number_I1_Resampled_10_Period,Bezugsfrequenz_I1_Resampled_10_Period,Signal_amp_I1_Resampled_10_Period] = Fast_Harmonics_1 (I1_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
[nyquist_number_I2_Resampled_10_Period,Bezugsfrequenz_I2_Resampled_10_Period,Signal_amp_I2_Resampled_10_Period] = Fast_Harmonics_1 (I2_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
[nyquist_number_I3_Resampled_10_Period,Bezugsfrequenz_I3_Resampled_10_Period,Signal_amp_I3_Resampled_10_Period] = Fast_Harmonics_1 (I3_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
[nyquist_number_U1_Resampled_10_Period,Bezugsfrequenz_U1_Resampled_10_Period,Signal_amp_U1_Resampled_10_Period] = Fast_Harmonics_1 (U1_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
[nyquist_number_U2_Resampled_10_Period,Bezugsfrequenz_U2_Resampled_10_Period,Signal_amp_U2_Resampled_10_Period] = Fast_Harmonics_1 (U2_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
[nyquist_number_U3_Resampled_10_Period,Bezugsfrequenz_U3_Resampled_10_Period,Signal_amp_U3_Resampled_10_Period] = Fast_Harmonics_1 (U3_Resample_10_Period_Final,F_Sample_10_Period_Resampled,T_S_10_Period_Resampled);
Figure_Show = 0;
if Figure_Show == 1
figure;
subplot(2,1,1);
plot([0:nyquist_number_data_I1]*Bezugsfrequenz_data_I1,Signal_amp_data_I1);ylabel('amplitude Currents [A]');title('Harmonic spectrom');hold on
plot([0:nyquist_number_data_I2]*Bezugsfrequenz_data_I2,Signal_amp_data_I2,'-.r');hold on
plot([0:nyquist_number_data_I3]*Bezugsfrequenz_data_I3,Signal_amp_data_I3,'--k');hold off
subplot(2,1,2);
plot([0:nyquist_number_data_U1]*Bezugsfrequenz_data_U1,Signal_amp_data_U1);xlabel('frequency [Hz]'),ylabel('amplitude Voltages [A]');hold on
plot([0:nyquist_number_data_U2]*Bezugsfrequenz_data_U2,Signal_amp_data_U2,'-.r');hold on
plot([0:nyquist_number_data_U3]*Bezugsfrequenz_data_U3,Signal_amp_data_U3,'.k','MarkerSize',0.1);hold off
end
% Grouping U1
data_specM_x_H = [0:nyquist_number_U1_Resampled_10_Period]*Bezugsfrequenz_U1_Resampled_10_Period;
deta_specM_x_V = data_specM_x_H';
data_specM_Y_V = Signal_amp_U1_Resampled_10_Period';
data_specM = data_specM_Y_V;
i = 11;
ii = 1;
while i < 502
harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
iTxt = (i-1)/10;
Results_U_H1(ii,jj) = harm_calc;
if i < 392
iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
Results_U_IH1(ii,jj) = iharm_calc;
end
if i <= 360
j1 = (i-1)*4+360+1;
j_end = j1+40;
hfharm_calc = 0;
while j1 <= j_end
hfharm_calc = hfharm_calc+power(data_specM(j1),2);
j1 =j1+1;
end
hfharm_calc = sqrt(hfharm_calc);
j1 = 200*(i-1)/10+1900;
Results_U_HF1(ii,jj) = hfharm_calc;
end
i = i+10;
ii = ii+1;
end
% Grouping U2
data_specM_x_H = [0:nyquist_number_U2_Resampled_10_Period]*Bezugsfrequenz_U2_Resampled_10_Period;
deta_specM_x_V = data_specM_x_H';
data_specM_Y_V = Signal_amp_U2_Resampled_10_Period';
data_specM = data_specM_Y_V;
i = 11;
ii = 1;
while i < 502
harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
iTxt = (i-1)/10;
Results_U_H2(ii,jj) = harm_calc;
if i < 392
iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
Results_U_IH2(ii,jj) = iharm_calc;
end
if i <= 360
j1 = (i-1)*4+360+1;
j_end = j1+40;
hfharm_calc = 0;
while j1 <= j_end
hfharm_calc = hfharm_calc+power(data_specM(j1),2);
j1 =j1+1;
end
hfharm_calc = sqrt(hfharm_calc);
j1 = 200*(i-1)/10+1900;
Results_U_HF2(ii,jj) = hfharm_calc;
end
i = i+10;
ii = ii+1;
end
% Grouping U3
data_specM_x_H = [0:nyquist_number_U3_Resampled_10_Period]*Bezugsfrequenz_U3_Resampled_10_Period;
deta_specM_x_V = data_specM_x_H';
data_specM_Y_V = Signal_amp_U3_Resampled_10_Period';
data_specM = data_specM_Y_V;
i = 11;
ii = 1;
while i < 502
harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
iTxt = (i-1)/10;
Results_U_H3(ii,jj) = harm_calc;
if i < 392
iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
Results_U_IH3(ii,jj) = iharm_calc;
end
if i <= 360
j1 = (i-1)*4+360+1;
j_end = j1+40;
hfharm_calc = 0;
while j1 <= j_end
hfharm_calc = hfharm_calc+power(data_specM(j1),2);
j1 =j1+1;
end
hfharm_calc = sqrt(hfharm_calc);
j1 = 200*(i-1)/10+1900;
Results_U_HF3(ii,jj) = hfharm_calc;
end
i = i+10;
ii = ii+1;
end
% Grouping I1
data_specM_x_H = [0:nyquist_number_I1_Resampled_10_Period]*Bezugsfrequenz_I1_Resampled_10_Period;
%deta_specM_x_V = data_specM_x_H';
data_specM_Y_V = Signal_amp_I1_Resampled_10_Period';
data_specM = data_specM_Y_V;
i = 11;
ii = 1;
while i < 502
harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
iTxt = (i-1)/10;
Results_I_H1(ii,jj) = harm_calc;
if i < 392
iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
Results_I_IH1(ii,jj) = iharm_calc;
end
if i <= 360
j1 = (i-1)*4+360+1;
j_end = j1+40;
hfharm_calc = 0;
while j1 <= j_end
hfharm_calc = hfharm_calc+power(data_specM(j1),2);
j1 =j1+1;
end
hfharm_calc = sqrt(hfharm_calc);
j1 = 200*(i-1)/10+1900;
Results_I_HF1(ii,jj) = hfharm_calc;
end
i = i+10;
ii = ii+1;
end
% Grouping I2
data_specM_x_H = [0:nyquist_number_I2_Resampled_10_Period]*Bezugsfrequenz_I2_Resampled_10_Period;
%deta_specM_x_V = data_specM_x_H';
data_specM_Y_V = Signal_amp_I2_Resampled_10_Period';
data_specM = data_specM_Y_V;
i = 11;
ii = 1;
while i < 502
harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
iTxt = (i-1)/10;
Results_I_H2(ii,jj) = harm_calc;
if i < 392
iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
Results_I_IH2(ii,jj) = iharm_calc;
end
if i <= 360
j1 = (i-1)*4+360+1;
j_end = j1+40;
hfharm_calc = 0;
while j1 <= j_end
hfharm_calc = hfharm_calc+power(data_specM(j1),2);
j1 =j1+1;
end
hfharm_calc = sqrt(hfharm_calc);
j1 = 200*(i-1)/10+1900;
Results_I_HF2(ii,jj) = hfharm_calc;
end
i = i+10;
ii = ii+1;
end
% Grouping I3
data_specM_x_H = [0:nyquist_number_I3_Resampled_10_Period]*Bezugsfrequenz_I3_Resampled_10_Period;
%deta_specM_x_V = data_specM_x_H';
data_specM_Y_V = Signal_amp_I3_Resampled_10_Period';
data_specM = data_specM_Y_V;
i = 11;
ii = 1;
while i < 502
harm_calc = sqrt(power(data_specM(i-1),2)+power(data_specM(i),2)+power(data_specM(i+1),2));
iTxt = (i-1)/10;
Results_I_H3(ii,jj) = harm_calc;
if i < 392
iharm_calc = sqrt(power(data_specM(i+2),2)+power(data_specM(i+3),2)+power(data_specM(i+4),2)+power(data_specM(i+5),2)+power(data_specM(i+6),2)+power(data_specM(i+7),2)+power(data_specM(i+8),2));
Results_I_IH3(ii,jj) = iharm_calc;
end
if i <= 360
j1 = (i-1)*4+360+1;
j_end = j1+40;
hfharm_calc = 0;
while j1 <= j_end
hfharm_calc = hfharm_calc+power(data_specM(j1),2);
j1 =j1+1;
end
hfharm_calc = sqrt(hfharm_calc);
j1 = 200*(i-1)/10+1900;
Results_I_HF3(ii,jj) = hfharm_calc;
end
i = i+10;
ii = ii+1;
end
if 0 == 1
v_alpha = sqrt(2)/3*(data_U1-0.5*data_U2-0.5*data_U3);
v_beta = sqrt(2)/3*(sqrt(3)/2*data_U2-sqrt(3)/2*data_U3);
I_alpha = sqrt(2)/3*(data_I1-0.5*data_I2-0.5*data_I3);
I_beta = sqrt(2)/3*(sqrt(3)/2*data_I2-sqrt(3)/2*data_I3);
P = 3*(v_alpha*I_alpha+v_beta*I_beta);
Q = 3*(-v_alpha*I_beta+v_beta*I_alpha);
Results_P(jj) = P;
Results_Q(jj) = Q;
end
j = j+20; % Calculation every 10 Periods
jj = jj+1;
end
%**********************************************************************
%----------------------------------END---------------------------------
%**********************************************************************
%-------------------------- Normalization -----------------------------
%**********************************************************************
Results_U_H1_pu = Results_U_H1/sqrt(2)/Un;
Results_U_H2_pu = Results_U_H2/sqrt(2)/Un;
Results_U_H3_pu = Results_U_H3/sqrt(2)/Un;
Results_I_H1_pu = Results_I_H1/sqrt(2)/In;
Results_I_H2_pu = Results_I_H2/sqrt(2)/In;
Results_I_H3_pu = Results_I_H3/sqrt(2)/In;
%----------------------------------------------------------------------
Results_U_IH1_pu = Results_U_IH1/sqrt(2)/Un;
Results_U_IH2_pu = Results_U_IH2/sqrt(2)/Un;
Results_U_IH3_pu = Results_U_IH3/sqrt(2)/Un;
Results_I_IH1_pu = Results_I_IH1/sqrt(2)/In;
Results_I_IH2_pu = Results_I_IH2/sqrt(2)/In;
Results_I_IH3_pu = Results_I_IH3/sqrt(2)/In;
%----------------------------------------------------------------------
Results_U_HF1_pu = Results_U_HF1/sqrt(2)/Un;
Results_U_HF2_pu = Results_U_HF2/sqrt(2)/Un;
Results_U_HF3_pu = Results_U_HF3/sqrt(2)/Un;
Results_I_HF1_pu = Results_I_HF1/sqrt(2)/In;
Results_I_HF2_pu = Results_I_HF2/sqrt(2)/In;
Results_I_HF3_pu = Results_I_HF3/sqrt(2)/In;
%**********************************************************************
%------------------------ End -----------------------------------------
%**********************************************************************
% ------------------- Mean and Max Value Calculation ------------------
% *********************************************************************
[~, colum] = size(Results_I_H1_pu);
if colum > 1
MeanHarm_U1 = mean(Results_U_H1_pu');
MeanHarm_U2 = mean(Results_U_H2_pu');
MeanHarm_U3 = mean(Results_U_H3_pu');
MeanHarm_I1 = mean(Results_I_H1_pu');
MeanHarm_I2 = mean(Results_I_H2_pu');
MeanHarm_I3 = mean(Results_I_H3_pu');
%----------------------------------------------------------------------
MeanIHarm_U1 = mean(Results_U_IH1_pu');
MeanIHarm_U2 = mean(Results_U_IH2_pu');
MeanIHarm_U3 = mean(Results_U_IH3_pu');
MeanIHarm_I1 = mean(Results_I_IH1_pu');
MeanIHarm_I2 = mean(Results_I_IH2_pu');
MeanIHarm_I3 = mean(Results_I_IH3_pu');
%----------------------------------------------------------------------
MeanHFHarm_U1 = mean(Results_U_HF1_pu');
MeanHFHarm_U2 = mean(Results_U_HF2_pu');
MeanHFHarm_U3 = mean(Results_U_HF3_pu');
MeanHFHarm_I1 = mean(Results_I_HF1_pu');
MeanHFHarm_I2 = mean(Results_I_HF2_pu');
MeanHFHarm_I3 = mean(Results_I_HF3_pu');
%----------------------------------------------------------------------
MaxHarm_U = max(max(MeanHarm_U1,MeanHarm_U2),MeanHarm_U3);
MaxHarm_I = max(max(MeanHarm_I1,MeanHarm_I2),MeanHarm_I3);
MaxIHarm_U = max(max(MeanIHarm_U1,MeanIHarm_U2),MeanIHarm_U3);
MaxIHarm_I = max(max(MeanIHarm_I1,MeanIHarm_I2),MeanIHarm_I3);
MaxHFHarm_U = max(max(MeanHFHarm_U1,MeanHFHarm_U2),MeanHFHarm_U3);
MaxHFHarm_I = max(max(MeanHFHarm_I1,MeanHFHarm_I2),MeanHFHarm_I3);
else
MeanHarm_I1 = Results_I_H1_pu;
MeanHarm_I2 = Results_I_H2_pu;
MeanHarm_I3 = Results_I_H3_pu;
MeanHarm_U1 = Results_U_H1_pu;
MeanHarm_U2 = Results_U_H2_pu;
MeanHarm_U3 = Results_U_H3_pu;
%----------------------------------------------------------------------
MeanIHarm_U1 = Results_U_IH1_pu;
MeanIHarm_U2 = Results_U_IH2_pu;
MeanIHarm_U3 = Results_U_IH3_pu;
MeanIHarm_I1 = Results_I_IH1_pu;
MeanIHarm_I2 = Results_I_IH2_pu;
MeanIHarm_I3 = Results_I_IH3_pu;
%----------------------------------------------------------------------
MeanHFHarm_U1 = Results_U_HF1_pu;
MeanHFHarm_U2 = Results_U_HF2_pu;
MeanHFHarm_U3 = Results_U_HF3_pu;
MeanHFHarm_I1 = Results_I_HF1_pu;
MeanHFHarm_I2 = Results_I_HF2_pu;
MeanHFHarm_I3 = Results_I_HF3_pu;
%----------------------------------------------------------------------
MaxHarm_U = max(max(MeanHarm_U1,MeanHarm_U2),MeanHarm_U3);
MaxHarm_I = max(max(MeanHarm_I1,MeanHarm_I2),MeanHarm_I3);
MaxIHarm_U = max(max(MeanIHarm_U1,MeanIHarm_U2),MeanIHarm_U3);
MaxIHarm_I = max(max(MeanIHarm_I1,MeanIHarm_I2),MeanIHarm_I3);
MaxHFHarm_U = max(max(MeanHFHarm_U1,MeanHFHarm_U2),MeanHFHarm_U3);
MaxHFHarm_I = max(max(MeanHFHarm_I1,MeanHFHarm_I2),MeanHFHarm_I3);
end
%**********************************************************************
%------------------------------- End ----------------------------------
%**********************************************************************
%--------------- Order (x-axie) Calculation and Ploting ---------------
%**********************************************************************
Order_Harm = 1:1:50;
Order_IHarm = 1:1:39;
%Order_HFHarm = 1:1:35;
%Order_HFHarm = 2100:1:2134;
Order_HFHarm = 2100:200:8900;
if figures
figure;
subplot(2,3,1)
bar(Order_Harm,MaxHarm_I);xlabel('Order Number');ylabel('Current Hamonics Levels I/In in pu');title('Current Harmonics spectrom')
subplot(2,3,2)
bar(Order_IHarm,MaxIHarm_I);xlabel('Order Number');ylabel('Current Hamonics Levels I/In in pu');title('Current interharmonics spectrom')
subplot(2,3,3)
bar(Order_HFHarm,MaxHFHarm_I);xlabel('Frequency in Hz');ylabel('Current Hamonics Levels I/In in pu');title('Current Higher Frequencies spectrom')
subplot(2,3,4)
bar(Order_Harm,MaxHarm_U);xlabel('Order Number');ylabel('Voltage Hamonics Levels U/Un in pu');title('Voltage Harmonics spectrom')
subplot(2,3,5)
bar(Order_IHarm,MaxIHarm_U);xlabel('Order Number');ylabel('Voltage Hamonics Levels U/Un in pu');title('Voltage interharmonics spectrom')
subplot(2,3,6)
bar(Order_HFHarm,MaxHFHarm_U);xlabel('Frequency in Hz');ylabel('Voltage Hamonics Levels U/Un in pu');title('Voltage Higher Frequencies spectrom')
end
%**********************************************************************
%------------------------------- End ----------------------------------
%**********************************************************************
%-------------------------- K, THD Factor  ----------------------------
%**********************************************************************
THD_U = sqrt(sum((MaxHarm_U(2:end).^2)'))/(MaxHarm_U(1));
THD_I = sqrt(sum((MaxHarm_I(2:end).^2)'))/(MaxHarm_I(1));
K_U = sqrt(sum((MaxHarm_U(2:end).^2)')/(MaxHarm_U(1)+sum((MaxHarm_U(2:end).^2)')));
K_I = sqrt(sum((MaxHarm_I(2:end).^2)')/(MaxHarm_I(1)+sum((MaxHarm_I(2:end).^2)')));
%**********************************************************************
%------------------------------- End ----------------------------------
%**********************************************************************
%------------------------ Saving Results ------------------------------
%**********************************************************************
save('Harmonics','THD_U','THD_I','K_U','K_I','MaxHarm_I','MaxHarm_U','MaxIHarm_I','MaxIHarm_U','MaxHFHarm_I','MaxHFHarm_U','MaxHarm_I','Order_IHarm','Order_Harm','Order_HFHarm')
%**********************************************************************
%------------------------------- End ----------------------------------
%**********************************************************************
end

```