IS_OCTAVE = exist(‘OCTAVE_VERSION’) ~= 0;

%% Check inputs

if (nargin ~= 3)

error(‘Invalid number of arguments’);

end

if (~isvector(u))

error(‘First input argument must be a vector’);

end

% convert to row vector if needed

u = reshape(u, 1, length(u));

if ((f_line ~= 50) && (f_line ~= 60))

error(‘Line frequency must be 50 or 60 Hz’);

end

if (fs < 2000)

warning(‘Sampling frequency should be >= 2000 Hz’);

end

%% Block 1: Input voltage adaptor

% remove DC component

u = u – mean(u);

% normalize input (scale with peak-amplitude value)

u_rms = sqrt(mean(u.^2));

u = u / (u_rms * sqrt(2));

%% Block 2: Quadratic demodulator

u_0 = u .^ 2;

%% Block 3: Bandpass and weighting filter

% bandpass filter

HIGHPASS_ORDER = 1;

HIGHPASS_CUTOFF = 0.05;

LOWPASS_ORDER = 6;

if (f_line == 50)

LOWPASS_CUTOFF = 35;

end

if (f_line == 60)

LOWPASS_CUTOFF = 42;

end

% subtract DC component to limit filter transients at start of simulation

u_0_ac = u_0 – mean(u_0);

[b_hp, a_hp] = butter(HIGHPASS_ORDER, HIGHPASS_CUTOFF / (fs / 2), ‘high’);

u_hp = filter(b_hp, a_hp, u_0_ac);

% smooth start of signal to avoid filter transient at start of simulation

smooth_limit = min(round(fs / 10), length(u_hp));

u_hp(1 : smooth_limit) = u_hp(1 : smooth_limit) .* linspace(0, 1, smooth_limit);

[b_bw, a_bw] = butter(LOWPASS_ORDER, LOWPASS_CUTOFF / (fs / 2), ‘low’);

u_bw = filter(b_bw, a_bw, u_hp);

% weighting filter

if (f_line == 50)

K = 1.74802;

LAMBDA = 2 * pi * 4.05981;

OMEGA1 = 2 * pi * 9.15494;

OMEGA2 = 2 * pi * 2.27979;

OMEGA3 = 2 * pi * 1.22535;

OMEGA4 = 2 * pi * 21.9;

end

if (f_line == 60)

K = 1.6357;

LAMBDA = 2 * pi * 4.167375;

OMEGA1 = 2 * pi * 9.077169;

OMEGA2 = 2 * pi * 2.939902;

OMEGA3 = 2 * pi * 1.394468;

OMEGA4 = 2 * pi * 17.31512;

end

num1 = ;

den1 = [1, 2 * LAMBDA, OMEGA1.^2];

num2 = ;

den2 = ;

if (IS_OCTAVE)

[b_w, a_w] = bilinear(conv(num1, num2), conv(den1, den2), 1 / fs);

else

[b_w, a_w] = bilinear(conv(num1, num2), conv(den1, den2), fs);

end

u_w = filter(b_w, a_w, u_bw);

%% Block 4: Squaring and smoothing

LOWPASS_2_ORDER = 1;

LOWPASS_2_CUTOFF = 1 / (2 * pi * 300e-3); % time constant 300 msec

SCALING_FACTOR = 1238400; % scaling of output to perceptibility scale (according [2])

u_q = u_w .^ 2;

[b_lp, a_lp] = butter(LOWPASS_2_ORDER, LOWPASS_2_CUTOFF / (fs / 2), ‘low’);

s = SCALING_FACTOR * filter(b_lp, a_lp, u_q);

%% Block 5: Statistical evaluation

NUMOF_CLASSES = 10000;

[bin_cnt, cpf.magnitude] = hist(s, NUMOF_CLASSES);

cpf.cum_probability = 100 * (1 – cumsum(bin_cnt) / sum(bin_cnt));

p_50s = mean([get_percentile(cpf, 30), get_percentile(cpf, 50), get_percentile(cpf, 80)]);

p_10s = mean([get_percentile(cpf, 6), get_percentile(cpf, 8), …

get_percentile(cpf, 10), get_percentile(cpf, 13), get_percentile(cpf, 17)]);

p_3s = mean([get_percentile(cpf, 2.2), get_percentile(cpf, 3), get_percentile(cpf, 4)]);

p_1s = mean([get_percentile(cpf, 0.7), get_percentile(cpf, 1), get_percentile(cpf, 1.5)]);

p_0_1 = get_percentile(cpf, 0.1);

P_st = sqrt(0.0314 * p_0_1 + 0.0525 * p_1s + 0.0657 * p_3s + …

۰٫۲۸ * p_10s + 0.08 * p_50s);

%% Optional graphical output

% time signals

if (SHOW_TIME_SIGNALS)

t = 0 : 1 / fs : (length(u) – 1) / fs;

filter_len = round(10 / 1000 * fs);

u_0_m = filter(ones(1, filter_len) / filter_len * 2, 1, u_0);

figure

clf

hold on

plot(t, u, ‘b’);

plot(t, u_0, ‘m’);

plot(t, u_hp, ‘r’);

plot(t, u_0_m, ‘c’);

hold off

legend(‘u’, ‘u_0’, ‘u_h_p’);

grid on

subplot(2, 2, 2)

hold on

plot(t, u_bw, ‘b’);

plot(t, u_w, ‘m’);

legend(‘u_b_w’, ‘u_w’);

hold off

grid on

subplot(2, 2, 3)

plot(t, u_q, ‘b’);

legend(‘u_q’);

grid on

subplot(2, 2, 4)