Warning: Signal Processing Toolbox is not installed. Skipping STFT analysis...
Example: MHKiT-MATLAB Power Module
This example will demonstrate the MHKiT Power Module functionality to compute power, instantaneous frequency, and harmonics from time series of voltage and current.
Load Power Data
We will begin by reading in time-series data of measured three phase (a,b,and,c) voltage and current. The IEC TS 62600-30 requires that you perform power quality assessments on a minimum of 10-min time-series data, but for this example we will only look at a fraction of that.
% Read in time-series data of voltage (V) and current (I)
power_table = readtable('./examples/data/power/2020224_181521_PowRaw.csv');
power_table
To use the MHKiT-MATLAB power module we need to create structures of current and voltage.
current.time = posixtime(power_table.Time_UTC); % setting time in current structure
voltage.time = posixtime(power_table.Time_UTC); % setting time in voltage structure
T2 = mergevars(power_table, [2 3 4]); % combining the voltage time series into one table variable
T3 = mergevars(T2, [3 4 5]); % combining the current time series into one table variable
current.current = T3.Var3;
voltage.voltage = T3.Var2;
Power Characteristics
The MHKiT power/characteristics submodule can be used to compute basic quantities of interest from voltage and current time series. In this example we will calculate active AC power and fundamental frequency, using the loaded voltage and current time-series.
To compute the active AC power, we will need the power factor. Power factor describes the efficiency of the energy device.
% Set the power factor for the system
power_factor = 0.96;
% Compute the instantaneous AC power in watts
ac_power = ac_power_three_phase(voltage, current, power_factor);
% Display the result
figure('Position', [100, 100, 1600, 600]);
plot(datetime(ac_power.time, "ConvertFrom", "posixtime"), ac_power.power);
title('AC Power');
ylabel('Power [W]');
xlabel('Time');
Fundamental Frequency
Using the 3 phase voltage measurements we can compute the fundamental frequency of the voltage time-series. The time-varying fundamental frequency is a required metric for power quality assessments. The function calc_fundamental_freq() provides two options of calculating the fundamental frequency: 1. Short-Time Fourier Transform (STFT) 2. Zero-Crossing Detection (ZCD) Here we illustrate the steps in both ways.
First, let's look at the STFT method. Note that, to obtain an accurate result, window settings for STFT is crucial, the user might need some exploratory analysis and try different sets of parameters.
Note: The STFT method requires the Signal Processing Toolbox
% Hilbert method: instantaneous_frequency
% inst_freq = instantaneous_frequency(voltage)
% STFT method:
% Check for Signal Processing Toolbox
if license('test', 'Signal_Toolbox')
% Set up method options
methodopts = {
'Window', int32(5000), ...
'OverlapLength', 3750, ...
'FFTLength', 50e3, ...
'FrequencyRange', 'onesided'
};
% Prep input and output
u_m = struct();
u_m.time = voltage.time;
fund_freq = struct();
fund_freq.time = voltage.time;
fund_freq.data = zeros(size(voltage.voltage));
for i = 1:3
u_m.data = voltage.voltage(:,i);
[~, freq] = calc_fundamental_freq(u_m, 'stft', methodopts);
fund_freq.data(:,i) = freq.data;
end
fund_freq
% Display the result
hold off;
figure('Position', [100, 100, 1600, 600]);
plot(datetime(fund_freq.time, "ConvertFrom", "posixtime"), fund_freq.data);
title('Fundamental Frequency: STFT');
ylabel('Frequency [Hz]');
xlabel('Time');
ylim([50, 70]);
else
warning('Signal Processing Toolbox is not installed. Skipping STFT analysis...');
end
Now, let's look at the ZCD method. Note that, ZCD method may not work for conditions of distorted voltage with multiple zero-crossings, like the one stated in IEC TS 61400-21-1 B.3.3. As stated above, exploratory analysis is needed.
% ZCD method:
% Prep input and output
u_m = struct();
u_m.time = voltage.time;
fund_freq = struct();
fund_freq.time = voltage.time;
fund_freq.data = zeros(size(voltage.voltage));
% For the zcd method, we don't need to worry about the method options
methodopts = {};
for i = 1:3
u_m.data = voltage.voltage(:,i);
[~, freq] = calc_fundamental_freq(u_m, 'zcd', methodopts);
fund_freq.data(:,i) = freq.data;
end
fund_freq
fund_freq = struct with fields:
time: [8000×1 double]
+.S14 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } 
+
Example: MHKiT-MATLAB Power Module
This example will demonstrate the MHKiT Power Module functionality to compute power, instantaneous frequency, and harmonics from time series of voltage and current.
Load Power Data
We will begin by reading in time-series data of measured three-phase (A, B, and C) voltage and current. IEC TS 62600-30 requires power quality assessments to use time series of at least 10 minutes, but for this example we will only look at a fraction of a second of data.
% Read in time-series data of voltage (V) and current (I)
power_table = readtable('./data/power/2020224_181521_PowRaw.csv');
% Display without scientific notation, then restore the default format
format shortG
power_table
format short
To use the MHKiT-MATLAB power module we need to create structures of current and voltage.
current.time = posixtime(power_table.Time_UTC); % setting time in current structure
voltage.time = posixtime(power_table.Time_UTC); % setting time in voltage structure
T2 = mergevars(power_table, [2 3 4]); % combining the voltage time series into one table variable
T3 = mergevars(T2, [3 4 5]); % combining the current time series into one table variable
current.current = T3.Var3;
voltage.voltage = T3.Var2;
IEC TS 62600-30 clause 6.1 requires the generator sign convention: power flowing from the marine energy converter to the grid is positive. In this dataset every sample of the instantaneous power (the sum of voltage times current over the three phases) is negative, which indicates the current transformers were connected with the opposite orientation. We therefore flip the sign of the measured currents so that the exported power is positive.
current.current = -current.current;
Before computing anything, let's look at the first 0.1 s of the three-phase voltage and current signals.
% Time since the start of the record [s]
elapsed_time = voltage.time - voltage.time(1);
first_tenth = elapsed_time <= 0.1;
figure('Position', [100, 100, 1600, 600]);
plot(elapsed_time(first_tenth), voltage.voltage(first_tenth, :));
title('Three-Phase Voltage');
xlabel('Time [s]');
ylabel('Voltage [V]');
ax = gca;
ax.YAxis.Exponent = 0; % show plain volts instead of scientific notation
legend('Phase A', 'Phase B', 'Phase C');
grid on;
figure('Position', [100, 100, 1600, 600]);
plot(elapsed_time(first_tenth), current.current(first_tenth, :));
title('Three-Phase Current');
xlabel('Time [s]');
ylabel('Current [A]');
legend('Phase A', 'Phase B', 'Phase C');
grid on;
Power Characteristics
The MHKiT power/characteristics submodule can be used to compute basic quantities of interest from voltage and current time series. In this example we will calculate active AC power and fundamental frequency, using the loaded voltage and current time-series.
To compute the active AC power, we will need the power factor, the ratio of active (real) power to apparent power (
, between 0 and 1). It describes how effectively the voltage and current at the connection deliver real power; it is not the device's conversion efficiency (
). For a grid-connected marine energy converter, the power factor at the grid connection is largely set by the power electronics' reactive power control (for example constant power factor, constant reactive power, or voltage control modes; IEC TS 62600-30 clause 6.7 and Table A.5). Here we assume a constant value of 0.96.
, between 0 and 1). It describes how effectively the voltage and current at the connection deliver real power; it is not the device's conversion efficiency (). For a grid-connected marine energy converter, the power factor at the grid connection is largely set by the power electronics' reactive power control (for example constant power factor, constant reactive power, or voltage control modes; IEC TS 62600-30 clause 6.7 and Table A.5). Here we assume a constant value of 0.96.% Set the power factor for the system
power_factor = 0.96;
% Compute the instantaneous AC power in watts
ac_power = ac_power_three_phase(voltage, current, power_factor);
% Display the result in kilowatts
figure('Position', [100, 100, 1600, 600]);
plot(datetime(ac_power.time, "ConvertFrom", "posixtime"), ac_power.power / 1000);
title('AC Power');
ylabel('Power [kW]');
xlabel('Time');
ax = gca;
ax.YAxis.Exponent = 0; % show plain kilowatts instead of scientific notation
Fundamental Frequency
Using the 3 phase voltage measurements we can compute the fundamental frequency of the voltage time series. The time-varying fundamental frequency is a required metric for power quality assessments. The function calc_fundamental_freq() provides two methods for calculating the fundamental frequency:
- Short-Time Fourier Transform (STFT)
- Zero-Crossing Detection (ZCD)
Here we illustrate both.
First, let's look at the STFT method. To obtain an accurate result the STFT window settings are crucial, so some exploratory analysis with different sets of parameters may be needed.
Note: The STFT method requires the Signal Processing Toolbox.
% Hilbert method: instantaneous_frequency
% inst_freq = instantaneous_frequency(voltage)
% STFT method:
% Check that the Signal Processing Toolbox functions are installed. license('test')
% only checks that a license exists, not that the toolbox is installed.
has_signal_toolbox = exist('stft', 'file') ~= 0 && exist('rectwin', 'file') ~= 0;
if has_signal_toolbox
% Set up method options
methodopts = {
'Window', int32(5000), ...
'OverlapLength', 3750, ...
'FFTLength', 50e3, ...
'FrequencyRange', 'onesided'
};
% Prep input and output
u_m = struct();
u_m.time = voltage.time;
fund_freq = struct();
fund_freq.time = voltage.time;
fund_freq.data = zeros(size(voltage.voltage));
for i = 1:3
u_m.data = voltage.voltage(:,i);
[~, freq] = calc_fundamental_freq(u_m, 'stft', methodopts);
fund_freq.data(:,i) = freq.data;
end
fund_freq
% Display the result
figure('Position', [100, 100, 1600, 600]);
plot(datetime(fund_freq.time, "ConvertFrom", "posixtime"), fund_freq.data);
title('Fundamental Frequency: STFT');
ylabel('Frequency [Hz]');
xlabel('Time');
ylim([50, 70]);
else
warning('Signal Processing Toolbox is not installed. Skipping STFT analysis...');
end
Now, let's look at the ZCD method. Note that the ZCD method may not work for distorted voltage with multiple zero crossings, such as the case described in IEC 61400-21-1:2019 Annex B.3.3. As stated above, exploratory analysis is needed.
% ZCD method:
% Prep input and output
u_m = struct();
u_m.time = voltage.time;
fund_freq = struct();
fund_freq.time = voltage.time;
fund_freq.data = zeros(size(voltage.voltage));
% For the zcd method, we don't need to worry about the method options
methodopts = {};
for i = 1:3
u_m.data = voltage.voltage(:,i);
[~, freq] = calc_fundamental_freq(u_m, 'zcd', methodopts);
fund_freq.data(:,i) = freq.data;
end
fund_freq
hold off;
figure('Position', [100, 100, 1600, 600]);
plot(datetime(fund_freq.time, "ConvertFrom", "posixtime"), fund_freq.data);
title('Fundamental Frequency: ZCD');
ylabel('Frequency [Hz]');
xlabel('Time');
ylim([50, 70]);
Power Quality
The power quality submodule can be used to compute voltage fluctuations, harmonics of current and voltage and current distortions following IEC TS 62600-30 and IEC 61000-4-7. Harmonics and harmonic distortion are required as part of a power quality assessment and characterize the stability of the power being produced. Specifically, we focus on the flicker coefficient for continuous operation (MV connected systems), which is a normalized measure of the flicker emission during continuous operation of the MEC unit.
% Set the sampling frequency of the dataset
sample_freq = 50000; % [Hz]
% Set the frequency of the grid the device would be connected to
grid_freq = 60; % [Hz]
% Set the rated current of the device
rated_current = 18.8; % [Amps]
% Calculate the current harmonics
h = harmonics(current, sample_freq, grid_freq);
% Display the results
figure('Position', [100, 100, 1600, 600]);
plot(h.harmonic, h.amplitude);
title('Current Harmonics');
xlabel('Frequency [Hz]');
ylabel('Harmonic Amplitude');
xlim([0, 900]);
Harmonic Subgroups
The harmonic subgroups calculations are based on IEC TS 62600-30. We can calculate them using our grid frequency and harmonics.
% Calculate Harmonic Subgroups
h_s = harmonic_subgroups(h, grid_freq);
Total Harmonic Current Distortion (THCD)
Compute the THCD from harmonic subgroups and rated current for the device.
% Finally we can compute the total harmonic current distortion as a percentage
% for each channel
THCD = total_harmonic_current_distortion(h_s, rated_current);
Flicker Assessment
Calculate the flicker coefficient following the steps in IEC TS 62600-30 section 7.2.1: MV connected marine energy converters.
1. Calculate the simulated voltage of the fictitious grid (u_fic). There are two options for calculation of u_fic.
First, the user can use flicker_ufic_workflow(), which wrapped up all steps needed to the calculation of u_fic from measured voltage (u_m) and current (i_m). All intermediate values can also be outputted for the user to debug.
Or, the user can call the corresponding functions sequentially. For illustration purpose, we show both ways.
Firstly, let's use flicker_ufic_workflow().
% Prep the input: look at the first phase for an example
u_m = struct();
u_m.time = voltage.time;
u_m.data = voltage.voltage(:,1);
i_m = struct();
i_m.time = current.time;
i_m.data = current.current(:,1);
Sr = 4.2e5; % rated power
Un = 14e3; % RMS value of the nominal voltage
SCR = 20; % short-circuit ratio
fg = 60; % nominal grid frequency
% Method settings for calculation of fundamental frequency
% of u_m, see above
method = 'zcd';
methodopts = {};
out = flicker_ufic_workflow(Sr, Un, SCR, fg, ...
u_m, i_m, method, methodopts);
u0 = out.u0;
out
figure('Position', [100, 100, 1600, 600]);
plot(datetime(fund_freq.time, "ConvertFrom", "posixtime"), fund_freq.data);
title('Fundamental Frequency: ZCD');
ylabel('Frequency [Hz]');
xlabel('Time');
ylim([50, 70]);
Power Quality
The power quality submodule can be used to compute voltage fluctuations (flicker), and harmonics and harmonic distortion of current and voltage, following IEC TS 62600-30 and IEC 61000-4-7. Harmonics and harmonic distortion are required as part of a power quality assessment and characterize the quality of the power being produced. We start with the current harmonics.
% Set the sampling frequency of the dataset
sample_freq = 50000; % [Hz]
% Set the frequency of the grid the device would be connected to
grid_freq = 60; % [Hz]
% Set the rated current of the device
rated_current = 18.8; % [A]
% Calculate the current harmonics
h = harmonics(current, sample_freq, grid_freq);
% Display the results
figure('Position', [100, 100, 1600, 600]);
plot(h.harmonic, h.amplitude);
title('Current Harmonics');
xlabel('Frequency [Hz]');
ylabel('Harmonic Amplitude [A]');
legend('Phase A', 'Phase B', 'Phase C');
xlim([0, 900]);
Harmonic Subgroups
IEC TS 62600-30 clause 7.3 requires the harmonic currents up to 50 times the grid frequency to be reported as harmonic subgroups. We calculate them from the harmonics and the grid frequency.
% Calculate harmonic subgroups
h_s = harmonic_subgroups(h, grid_freq);
% Display the first ten harmonic orders (row 1 is DC, row 2 is the fundamental)
harmonic_order = (0:9)';
array2table([harmonic_order, h_s.harmonic(1:10), h_s.amplitude(1:10, :)], ...
'VariableNames', {'Order', 'Frequency [Hz]', 'Phase A [A]', 'Phase B [A]', 'Phase C [A]'})
Total Harmonic Current Distortion (THCD)
Finally we compute the total harmonic current distortion as a percentage. Matching MHKiT-Python, the subgrouped harmonics of order 2 to 49 are summed and normalized by the fundamental subgroup, and the result is for the first column (Phase A) only. The rated current is accepted but not yet used; the IEC TS 62600-30 formula (7) normalization by rated current and per-phase output are under review.
THCD = total_harmonic_current_distortion(h_s, rated_current)
Flicker Assessment
Calculate the flicker coefficient following the steps in IEC TS 62600-30 section 7.2.1: MV connected marine energy converters.
The first step is to calculate the simulated voltage of the fictitious grid (u_fic). There are two options for calculating u_fic.
First, the user can use flicker_ufic_workflow(), which wraps up all the steps needed to calculate u_fic from the measured voltage (u_m) and current (i_m). All intermediate values are also returned to help with debugging.
Alternatively, the user can call the corresponding functions sequentially. For illustration, we show both ways.
Firstly, let's use flicker_ufic_workflow().
% Prep the input: look at the first phase for an example
u_m = struct();
u_m.time = voltage.time;
u_m.data = voltage.voltage(:,1);
i_m = struct();
i_m.time = current.time;
i_m.data = current.current(:,1);
Sr = 4.2e5; % rated apparent power [VA]
Un = 14e3; % RMS value of the nominal voltage [V]
SCR = 20; % short-circuit ratio
fg = 60; % nominal grid frequency [Hz]
% Method settings for calculation of fundamental frequency
% of u_m, see above
method = 'zcd';
methodopts = {};
out = flicker_ufic_workflow(Sr, Un, SCR, fg, ...
u_m, i_m, method, methodopts);
u0 = out.u0;
out
Calculating u_fic can be tricky due to the part of calculating the ideal voltage source (u0), as mentioned in the IEC TS 62600-30, u0 should fulfill the following two requirements:
(1) be without any fluctuations (2) have the same electrical angle as the fundamental of u_m
Therefore, here we compare u0 and u_m to check if our output fulfills these requirements.
% Check u0
hold off;
figure('Position', [100, 100, 1600, 600]);
plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u_m.data);
hold on;
figure('Position', [100, 100, 1600, 600]);
plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u0);
legend('u_m', 'u_0');
Secondly, we can also calculate u_fic step by step. The calculation of u_fic can be achieved sequentially in 5 steps:
Step 1. Construct the fictitious grid: calculate resistance and inductance using calc_Rfic_Lfic().
[Rfic, Lfic] = calc_Rfic_Lfic(Sr, SCR, Un, fg);
Step 2. Calculate the fundamental frequency and alpha_0: -opt1: method ZCD: method = 'ZCD'; methodopts={} -opt2: method STFT: method = 'stft'; methodopts = {'Window',rectwin(M),... 'OverlapLength',L,'FFTLength',128,'FrequencyRange','onesided'}
if license('test', 'Signal_Toolbox')
method = 'stft';
methodopts = {
'Window', rectwin(5000), ...
'OverlapLength', 3750, ...
'FFTLength', 50e3, ...
'FrequencyRange', 'onesided'
};
else
warning('Signal Processing Toolbox is not installed. Using ZCD method...');
method = 'ZCD';
methodopts = {};
end
[alpha0, freq] = calc_fundamental_freq(u_m, method, methodopts);
Step 3. Calculate the electrical angle (alpha_m) of the fundamental of u_m using calc_electrical_angle().
alpha_m = calc_electrical_angle(freq, alpha0);
Step 4. Calculate u0 from alpha_m and nominal voltage (Un) using calc_ideal_voltage().
u0 = calc_ideal_voltage(Un, alpha_m);
Step 5. Calculate u_fic using calc_simulated_voltage().
u_fic = calc_simulated_voltage(u0, i_m, Rfic, Lfic);
disp(size(u_fic));
Again, we check our derived u0 to see if it fulfills the requirements.
% Check u0
hold off;
figure('Position', [100, 100, 1600, 600]);
plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u_m.data);
hold on;
figure('Position', [100, 100, 1600, 600]);
plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u0);
legend('u_m', 'u_0');
Calculate the flicker emission value (P_stfic) from u_fic
Input u_fic into an appropriate digital flickermeter to get the P_stfic.
According to the standard, the evaluation time is 10 min, therefore, to fulfill this requirement, we need data length >= 10min. And note that we may need to discard the first 20s depending on the performance of the flickermeter we choose.
% Prep input for the digital flickermeter
time = u_m.time - u_m.time(1);
u_fic_in = timeseries(u_fic(:,2), time);
u_fic_in
Calculating u_fic can be tricky because of the ideal voltage source (u0). As stated in IEC TS 62600-30 clause 7.2.2, u0 should fulfill the following two requirements:
(1) be without any fluctuations (2) have the same electrical angle as the fundamental of u_m
Therefore, we compare u0 and u_m to check that our output fulfills these requirements.
% Check u0
figure('Position', [100, 100, 1600, 600]);
plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u_m.data);
hold on;
plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u0);
hold off;
legend('u_m', 'u_0');
xlabel('Time');
ylabel('Voltage [V]');
Secondly, we can also calculate u_fic step by step. The calculation of u_fic can be achieved sequentially in 5 steps:
Step 1. Construct the fictitious grid: calculate resistance and inductance using calc_Rfic_Lfic().
[Rfic, Lfic] = calc_Rfic_Lfic(Sr, SCR, Un, fg);
Step 2. Calculate the fundamental frequency and alpha_0: -opt1: method ZCD: method = 'ZCD'; methodopts={} -opt2: method STFT: method = 'stft'; methodopts = {'Window',rectwin(M),... 'OverlapLength',L,'FFTLength',128,'FrequencyRange','onesided'}
if has_signal_toolbox
method = 'stft';
methodopts = {
'Window', rectwin(5000), ...
'OverlapLength', 3750, ...
'FFTLength', 50e3, ...
'FrequencyRange', 'onesided'
};
else
warning('Signal Processing Toolbox is not installed. Using ZCD method...');
method = 'ZCD';
methodopts = {};
end
[alpha0, freq] = calc_fundamental_freq(u_m, method, methodopts);
Step 3. Calculate the electrical angle (alpha_m) of the fundamental of u_m using calc_electrical_angle().
alpha_m = calc_electrical_angle(freq, alpha0);
Step 4. Calculate u0 from alpha_m and nominal voltage (Un) using calc_ideal_voltage().
u0 = calc_ideal_voltage(Un, alpha_m);
Step 5. Calculate u_fic using calc_simulated_voltage().
u_fic = calc_simulated_voltage(u0, i_m, Rfic, Lfic);
Again, we check our derived u0 to see if it fulfills the requirements.
% Check u0
figure('Position', [100, 100, 1600, 600]);
plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u_m.data);
hold on;
plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u0);
hold off;
legend('u_m', 'u_0');
xlabel('Time');
ylabel('Voltage [V]');
Calculate the flicker emission value (P_stfic) from u_fic
Input u_fic into an appropriate digital flickermeter to get P_stfic.
According to the standard, the evaluation time is 10 min, so the data must be at least 10 min long. Depending on the performance of the chosen flickermeter, the first 20 s may need to be discarded.
u_fic has one column for each fictitious grid impedance phase angle, psi_k = 30, 50, 70, and 85 degrees (see calc_Rfic_Lfic()). IEC TS 62600-30 requires the flicker coefficient to be reported for all four angles (Annex A, Table A.6); here we illustrate a single angle, psi_k = 50 degrees.
% Select the fictitious grid impedance phase angle
psi_k = [30, 50, 70, 85]; % [degrees], column order of u_fic, Rfic, and Lfic
psi_idx = find(psi_k == 50);
% Prep input for the digital flickermeter
time = u_m.time - u_m.time(1);
u_fic_in = timeseries(u_fic(:, psi_idx), time);
u_fic_in
A digital flickermeter that also contains an interface of statistical analysis in simulink can be used to generate the instantaneous flicker level (out.S5) and then calculate the Pst using the statistical tool associated with it. More details can be found here.
Note that the user needs to specify RMS value for voltage, sampling frequencies, and other parameters to achieve a valid result. In addition, there may be some unreliable peaks in the instantaneous flicker levels derived from this particular flickermeter, when performing statistical analysis, be sure to discard the S5 data at the very beginning of it.
% Prep input for the statistical tool
% out.S5 is the result generated from the digital flickermeter.
% S5 = out.S5;
After getting the instantaneous flicker levels, double-click the digital flickermeter to open statistical tools to calculate P_stfic.
Calculate the flicker coefficient
Provide a flickermeter output (P_stfic). In this example, data duration is < 1s, so we fake a number here for P_stfic
P_stfic = 0.10016;
Xfic = 2*pi*fg*Lfic;
S_kfic = (Un^2)./sqrt(Rfic.^2 + Xfic.^2);
coef_flicker = calc_flicker_coefficient(P_stfic, S_kfic(1), Sr);
coef_flicker
coef_flicker = 2.0032
Test the performance of the flicker assessment process
This part of the example illustrates how to test the performance of the user's flicker assessment process, including methods to generate simulated voltage of the fictitious grid and the performance of the (digital) flickermeter, following the guidance provided in the IEC TS 61400-21-1 Annex B.3: Verification test of the measurement procedure for flicker.
The verification is straightforward, we generate testing data with predetermined flicker coefficients and verify the flicker assessment process by comparing the derived flicker coefficient with the predetermined ones. The testing data generated will be measured voltage and current time series (u_m and i_m), and the user can use gen_test_data() to generate test data for five different scenarios: (1) pure sine wave (coef=0) and scenarios (2)-(5) as stated Annex B.3.2 - B.3.5, and test them one by one. Here we illustrate the usage of gen_test_data by generating distorted u_m(t) with multiple zero crossings (see Annex B.3.3 for details).
% Set up proper input parameters:
% scenario 3 with duration of 10s
opt = 3; T=10;
% rated parameters:
Un=12e3; In=144; fs=50e3; Sr = 3e6;
% other parameters for generating distorted data
fg=60; SCR=20; fm=25;
fv=0.5; % fv only needed for scenario B.3.4
% According to IEC TS 61400-21-1 Table B.2:
DeltaI_I = [4.763 5.726 7.640 9.488];
[i_m,u_m]=gen_test_data(Un,In,fg,fs,fm,fv,DeltaI_I,opt,T)
A digital flickermeter that also contains an interface of statistical analysis in simulink can be used to generate the instantaneous flicker level (out.S5) and then calculate the Pst using the statistical tool associated with it. More details can be found here.
Note that the user needs to specify RMS value for voltage, sampling frequencies, and other parameters to achieve a valid result. In addition, there may be some unreliable peaks in the instantaneous flicker levels derived from this particular flickermeter, when performing statistical analysis, be sure to discard the S5 data at the very beginning of it.
% Prep input for the statistical tool
% out.S5 is the result generated from the digital flickermeter.
% S5 = out.S5;
After getting the instantaneous flicker levels, double-click the digital flickermeter to open statistical tools to calculate P_stfic.
Calculate the flicker coefficient
Provide the flickermeter output (P_stfic) for the selected impedance phase angle. In this example the data duration is less than 1 s, so we use a placeholder value for P_stfic.
P_stfic = 0.10016;
Xfic = 2*pi*fg*Lfic;
S_kfic = (Un^2)./sqrt(Rfic.^2 + Xfic.^2); % IEC TS 62600-30 formula (5)
coef_flicker = calc_flicker_coefficient(P_stfic, S_kfic(psi_idx), Sr)
Test the performance of the flicker assessment process
This part of the example illustrates how to test the performance of the user's flicker assessment process, including methods to generate simulated voltage of the fictitious grid and the performance of the (digital) flickermeter, following the guidance provided in IEC 61400-21-1:2019 Annex B.3: Verification test of the measurement procedure for flicker.
The verification is straightforward: we generate test data with predetermined flicker coefficients and verify the flicker assessment process by comparing the derived flicker coefficients with the predetermined ones. The generated test data are measured voltage and current time series (u_m and i_m). The user can use gen_test_data() to generate test data for five scenarios: (1) a pure sine wave (coefficient = 0) and (2)-(5) the scenarios in Annex B.3.2 to B.3.5, and test them one by one. Here we illustrate gen_test_data() by generating a distorted u_m(t) with multiple zero crossings (Annex B.3.3).
% Set up proper input parameters:
% scenario 3 with duration of 10s
opt = 3; T=10;
% rated parameters: voltage [V], current [A], sample rate [Hz], apparent power [VA]
Un=12e3; In=144; fs=50e3; Sr = 3e6;
% other parameters for generating distorted data
fg=60; SCR=20; fm=25;
fv=0.5; % fv only needed for scenario B.3.4
% According to IEC 61400-21-1:2019 Table B.2:
DeltaI_I = [4.763 5.726 7.640 9.488];
[i_m,u_m]=gen_test_data(Un,In,fg,fs,fm,fv,DeltaI_I,opt,T)
The generated i_m(ntime, nphi_k) has 4 time-series, corresponding to phi_k = 30, 50, 70, 85. While u_m is always the same, thus has only one time series.
figure('Position', [100, 100, 1600, 600]);
hold off;
plot(i_m.time, i_m.data(:,1));
xlim([0 1]);
xlabel("Time");
ylabel("i_m");
figure('Position', [100, 100, 1600, 600]);
plot(u_m.time, u_m.data);
xlim([0 0.1]);
xlabel("Time");
ylabel("u_m");
grid on;
The generated i_m has four time series (columns), one for each impedance phase angle psi_k = 30, 50, 70, and 85 degrees. u_m is the same for all four angles, so it has a single time series.
figure('Position', [100, 100, 1600, 600]);
plot(i_m.time, i_m.data(:,1));
xlim([0 1]);
xlabel("Time [s]");
ylabel("i_m [A]");
title("Generated current, \psi_k = 30 degrees");
figure('Position', [100, 100, 1600, 600]);
plot(u_m.time, u_m.data);
xlim([0 0.1]);
xlabel("Time [s]");
ylabel("u_m [V]");
title("Generated voltage with multiple zero crossings");
grid on;
diff --git a/examples/power_example.m b/examples/power_example.m index 8d566e8ae..d4b807865 100644 --- a/examples/power_example.m +++ b/examples/power_example.m @@ -5,14 +5,18 @@ % and harmonics from time series of voltage and current. %% Load Power Data -% We will begin by reading in time-series data of measured three phase (a,b,and,c) -% voltage and current. The IEC TS 62600-30 requires that you perform power quality -% assessments on a minimum of 10-min time-series data, but for this example we -% will only look at a fraction of that. +% We will begin by reading in time-series data of measured three-phase (A, B, +% and C) voltage and current. IEC TS 62600-30 requires power quality assessments +% to use time series of at least 10 minutes, but for this example we will only +% look at a fraction of a second of data. % Read in time-series data of voltage (V) and current (I) -power_table = readtable('./examples/data/power/2020224_181521_PowRaw.csv'); +power_table = readtable('./data/power/2020224_181521_PowRaw.csv'); + +% Display without scientific notation, then restore the default format +format shortG power_table +format short %% % To use the MHKiT-MATLAB power module we need to create structures of current @@ -25,6 +29,42 @@ current.current = T3.Var3; voltage.voltage = T3.Var2; +%% +% IEC TS 62600-30 clause 6.1 requires the generator sign convention: power +% flowing from the marine energy converter to the grid is positive. In this +% dataset every sample of the instantaneous power (the sum of voltage times +% current over the three phases) is negative, which indicates the current +% transformers were connected with the opposite orientation. We therefore +% flip the sign of the measured currents so that the exported power is positive. + +current.current = -current.current; + +%% +% Before computing anything, let's look at the first 0.1 s of the three-phase +% voltage and current signals. + +% Time since the start of the record [s] +elapsed_time = voltage.time - voltage.time(1); +first_tenth = elapsed_time <= 0.1; + +figure('Position', [100, 100, 1600, 600]); +plot(elapsed_time(first_tenth), voltage.voltage(first_tenth, :)); +title('Three-Phase Voltage'); +xlabel('Time [s]'); +ylabel('Voltage [V]'); +ax = gca; +ax.YAxis.Exponent = 0; % show plain volts instead of scientific notation +legend('Phase A', 'Phase B', 'Phase C'); +grid on; + +figure('Position', [100, 100, 1600, 600]); +plot(elapsed_time(first_tenth), current.current(first_tenth, :)); +title('Three-Phase Current'); +xlabel('Time [s]'); +ylabel('Current [A]'); +legend('Phase A', 'Phase B', 'Phase C'); +grid on; + %% Power Characteristics % % The MHKiT |power/characteristics| submodule can be used to compute basic quantities @@ -32,8 +72,15 @@ % active AC power and fundamental frequency, using the loaded voltage and current % time-series. % -% To compute the active AC power, we will need the power factor. Power factor -% describes the efficiency of the energy device. +% To compute the active AC power, we will need the power factor, the ratio of +% active (real) power to apparent power ($PF = P/S$, between 0 and 1). It +% describes how effectively the voltage and current at the connection deliver +% real power; it is not the device's conversion efficiency ($P_{out}/P_{in}$). +% For a grid-connected marine energy converter, the power factor at the grid +% connection is largely set by the power electronics' reactive power control +% (for example constant power factor, constant reactive power, or voltage +% control modes; IEC TS 62600-30 clause 6.7 and Table A.5). Here we assume a +% constant value of 0.96. % Set the power factor for the system power_factor = 0.96; @@ -41,36 +88,42 @@ % Compute the instantaneous AC power in watts ac_power = ac_power_three_phase(voltage, current, power_factor); -% Display the result +% Display the result in kilowatts figure('Position', [100, 100, 1600, 600]); -plot(datetime(ac_power.time, "ConvertFrom", "posixtime"), ac_power.power); +plot(datetime(ac_power.time, "ConvertFrom", "posixtime"), ac_power.power / 1000); title('AC Power'); -ylabel('Power [W]'); +ylabel('Power [kW]'); xlabel('Time'); +ax = gca; +ax.YAxis.Exponent = 0; % show plain kilowatts instead of scientific notation %% % *Fundamental Frequency* % % Using the 3 phase voltage measurements we can compute the fundamental frequency -% of the voltage time-series. The time-varying fundamental frequency is a required +% of the voltage time series. The time-varying fundamental frequency is a required % metric for power quality assessments. The function |calc_fundamental_freq()| -% provides two options of calculating the fundamental frequency: -% 1. Short-Time Fourier Transform (STFT) -% 2. Zero-Crossing Detection (ZCD) -% Here we illustrate the steps in both ways. +% provides two methods for calculating the fundamental frequency: +% +% # Short-Time Fourier Transform (STFT) +% # Zero-Crossing Detection (ZCD) +% +% Here we illustrate both. % -% First, let's look at the STFT method. Note that, to obtain an accurate result, -% window settings for STFT is crucial, the user might need some exploratory analysis -% and try different sets of parameters. +% First, let's look at the STFT method. To obtain an accurate result the STFT +% window settings are crucial, so some exploratory analysis with different sets +% of parameters may be needed. % -% Note: The STFT method requires the Signal Processing Toolbox +% Note: The STFT method requires the Signal Processing Toolbox. % Hilbert method: instantaneous_frequency % inst_freq = instantaneous_frequency(voltage) % STFT method: -% Check for Signal Processing Toolbox -if license('test', 'Signal_Toolbox') +% Check that the Signal Processing Toolbox functions are installed. license('test') +% only checks that a license exists, not that the toolbox is installed. +has_signal_toolbox = exist('stft', 'file') ~= 0 && exist('rectwin', 'file') ~= 0; +if has_signal_toolbox % Set up method options methodopts = { 'Window', int32(5000), ... @@ -94,7 +147,6 @@ fund_freq % Display the result - hold off; figure('Position', [100, 100, 1600, 600]); plot(datetime(fund_freq.time, "ConvertFrom", "posixtime"), fund_freq.data); title('Fundamental Frequency: STFT'); @@ -106,9 +158,9 @@ end %% -% Now, let's look at the ZCD method. Note that, ZCD method may not work for -% conditions of distorted voltage with multiple zero-crossings, like the one stated -% in IEC TS 61400-21-1 B.3.3. As stated above, exploratory analysis is needed. +% Now, let's look at the ZCD method. Note that the ZCD method may not work for +% distorted voltage with multiple zero crossings, such as the case described in +% IEC 61400-21-1:2019 Annex B.3.3. As stated above, exploratory analysis is needed. % ZCD method: % Prep input and output @@ -128,8 +180,6 @@ end fund_freq -hold off; - figure('Position', [100, 100, 1600, 600]); plot(datetime(fund_freq.time, "ConvertFrom", "posixtime"), fund_freq.data); title('Fundamental Frequency: ZCD'); @@ -139,13 +189,11 @@ %% Power Quality % -% The power quality submodule can be used to compute voltage fluctuations, harmonics -% of current and voltage and current distortions following IEC TS 62600-30 and -% IEC 61000-4-7. Harmonics and harmonic distortion are required as part of a power -% quality assessment and characterize the stability of the power being produced. -% Specifically, we focus on the flicker coefficient for continuous operation (MV -% connected systems), which is a normalized measure of the flicker emission during -% continuous operation of the MEC unit. +% The power quality submodule can be used to compute voltage fluctuations +% (flicker), and harmonics and harmonic distortion of current and voltage, +% following IEC TS 62600-30 and IEC 61000-4-7. Harmonics and harmonic distortion +% are required as part of a power quality assessment and characterize the +% quality of the power being produced. We start with the current harmonics. % Set the sampling frequency of the dataset sample_freq = 50000; % [Hz] @@ -154,7 +202,7 @@ grid_freq = 60; % [Hz] % Set the rated current of the device -rated_current = 18.8; % [Amps] +rated_current = 18.8; % [A] % Calculate the current harmonics h = harmonics(current, sample_freq, grid_freq); @@ -164,27 +212,36 @@ plot(h.harmonic, h.amplitude); title('Current Harmonics'); xlabel('Frequency [Hz]'); -ylabel('Harmonic Amplitude'); +ylabel('Harmonic Amplitude [A]'); +legend('Phase A', 'Phase B', 'Phase C'); xlim([0, 900]); %% % *Harmonic Subgroups* % -% The harmonic subgroups calculations are based on IEC TS 62600-30. We can calculate -% them using our grid frequency and harmonics. +% IEC TS 62600-30 clause 7.3 requires the harmonic currents up to 50 times the +% grid frequency to be reported as harmonic subgroups. We calculate them from +% the harmonics and the grid frequency. -% Calculate Harmonic Subgroups +% Calculate harmonic subgroups h_s = harmonic_subgroups(h, grid_freq); +% Display the first ten harmonic orders (row 1 is DC, row 2 is the fundamental) +harmonic_order = (0:9)'; +array2table([harmonic_order, h_s.harmonic(1:10), h_s.amplitude(1:10, :)], ... + 'VariableNames', {'Order', 'Frequency [Hz]', 'Phase A [A]', 'Phase B [A]', 'Phase C [A]'}) + %% % *Total Harmonic Current Distortion (THCD)* % -% Compute the THCD from harmonic subgroups and rated current for the device. - -% Finally we can compute the total harmonic current distortion as a percentage -% for each channel +% Finally we compute the total harmonic current distortion as a percentage. +% Matching MHKiT-Python, the subgrouped harmonics of order 2 to 49 are summed +% and normalized by the fundamental subgroup, and the result is for the first +% column (Phase A) only. The rated current is accepted but not yet used; the +% IEC TS 62600-30 formula (7) normalization by rated current and per-phase +% output are under review. -THCD = total_harmonic_current_distortion(h_s, rated_current); +THCD = total_harmonic_current_distortion(h_s, rated_current) %% % *Flicker Assessment* @@ -192,15 +249,15 @@ % Calculate the flicker coefficient following the steps in IEC TS 62600-30 section % 7.2.1: MV connected marine energy converters. % -% 1. Calculate the simulated voltage of the fictitious grid (u_fic). -% There are two options for calculation of u_fic. +% The first step is to calculate the simulated voltage of the fictitious grid +% (u_fic). There are two options for calculating u_fic. % -% First, the user can use |flicker_ufic_workflow()|, which wrapped up all steps -% needed to the calculation of u_fic from measured voltage (u_m) and current (i_m). -% All intermediate values can also be outputted for the user to debug. +% First, the user can use |flicker_ufic_workflow()|, which wraps up all the steps +% needed to calculate u_fic from the measured voltage (u_m) and current (i_m). +% All intermediate values are also returned to help with debugging. % -% Or, the user can call the corresponding functions sequentially. For illustration -% purpose, we show both ways. +% Alternatively, the user can call the corresponding functions sequentially. For +% illustration, we show both ways. % % Firstly, let's use |flicker_ufic_workflow()|. @@ -212,10 +269,10 @@ i_m.time = current.time; i_m.data = current.current(:,1); -Sr = 4.2e5; % rated power -Un = 14e3; % RMS value of the nominal voltage +Sr = 4.2e5; % rated apparent power [VA] +Un = 14e3; % RMS value of the nominal voltage [V] SCR = 20; % short-circuit ratio -fg = 60; % nominal grid frequency +fg = 60; % nominal grid frequency [Hz] % Method settings for calculation of fundamental frequency % of u_m, see above @@ -228,24 +285,25 @@ out %% -% Calculating u_fic can be tricky due to the part of calculating the ideal voltage -% source (u0), as mentioned in the IEC TS 62600-30, u0 should fulfill the following -% two requirements: +% Calculating u_fic can be tricky because of the ideal voltage source (u0). As +% stated in IEC TS 62600-30 clause 7.2.2, u0 should fulfill the following two +% requirements: % % (1) be without any fluctuations % (2) have the same electrical angle as the fundamental of u_m % -% Therefore, here we compare u0 and u_m to check if our output fulfills these +% Therefore, we compare u0 and u_m to check that our output fulfills these % requirements. % Check u0 -hold off; figure('Position', [100, 100, 1600, 600]); plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u_m.data); hold on; -figure('Position', [100, 100, 1600, 600]); plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u0); +hold off; legend('u_m', 'u_0'); +xlabel('Time'); +ylabel('Voltage [V]'); %% % Secondly, we can also calculate u_fic step by step. The calculation of u_fic @@ -263,7 +321,7 @@ % -opt2: method STFT: method = 'stft'; methodopts = {'Window',rectwin(M),... % 'OverlapLength',L,'FFTLength',128,'FrequencyRange','onesided'} -if license('test', 'Signal_Toolbox') +if has_signal_toolbox method = 'stft'; methodopts = { 'Window', rectwin(5000), ... @@ -294,32 +352,41 @@ % Step 5. Calculate u_fic using |calc_simulated_voltage()|. u_fic = calc_simulated_voltage(u0, i_m, Rfic, Lfic); -disp(size(u_fic)); %% % Again, we check our derived u0 to see if it fulfills the requirements. % Check u0 -hold off; figure('Position', [100, 100, 1600, 600]); plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u_m.data); hold on; -figure('Position', [100, 100, 1600, 600]); plot(datetime(u_m.time, "ConvertFrom", "posixtime"), u0); +hold off; legend('u_m', 'u_0'); +xlabel('Time'); +ylabel('Voltage [V]'); %% % *Calculate the flicker emission value (P_stfic) from u_fic* % -% Input u_fic into an appropriate digital flickermeter to get the P_stfic. +% Input u_fic into an appropriate digital flickermeter to get P_stfic. % -% According to the standard, the evaluation time is 10 min, therefore, to fulfill -% this requirement, we need data length >= 10min. And note that we may need to -% discard the first 20s depending on the performance of the flickermeter we choose. +% According to the standard, the evaluation time is 10 min, so the data must be +% at least 10 min long. Depending on the performance of the chosen flickermeter, +% the first 20 s may need to be discarded. +% +% u_fic has one column for each fictitious grid impedance phase angle, psi_k = +% 30, 50, 70, and 85 degrees (see |calc_Rfic_Lfic()|). IEC TS 62600-30 requires +% the flicker coefficient to be reported for all four angles (Annex A, Table +% A.6); here we illustrate a single angle, psi_k = 50 degrees. + +% Select the fictitious grid impedance phase angle +psi_k = [30, 50, 70, 85]; % [degrees], column order of u_fic, Rfic, and Lfic +psi_idx = find(psi_k == 50); % Prep input for the digital flickermeter time = u_m.time - u_m.time(1); -u_fic_in = timeseries(u_fic(:,2), time); +u_fic_in = timeseries(u_fic(:, psi_idx), time); u_fic_in %% @@ -346,15 +413,14 @@ %% % *Calculate the flicker coefficient* % -% Provide a flickermeter output (P_stfic). -% In this example, data duration is < 1s, -% so we fake a number here for P_stfic +% Provide the flickermeter output (P_stfic) for the selected impedance phase +% angle. In this example the data duration is less than 1 s, so we use a +% placeholder value for P_stfic. P_stfic = 0.10016; Xfic = 2*pi*fg*Lfic; -S_kfic = (Un^2)./sqrt(Rfic.^2 + Xfic.^2); -coef_flicker = calc_flicker_coefficient(P_stfic, S_kfic(1), Sr); -coef_flicker +S_kfic = (Un^2)./sqrt(Rfic.^2 + Xfic.^2); % IEC TS 62600-30 formula (5) +coef_flicker = calc_flicker_coefficient(P_stfic, S_kfic(psi_idx), Sr) %% % *Test the performance of the flicker assessment process* @@ -362,45 +428,46 @@ % This part of the example illustrates how to test the performance of the user's % flicker assessment process, including methods to generate simulated voltage % of the fictitious grid and the performance of the (digital) flickermeter, following -% the guidance provided in the IEC TS 61400-21-1 Annex B.3: Verification test of +% the guidance provided in IEC 61400-21-1:2019 Annex B.3: Verification test of % the measurement procedure for flicker. % -% The verification is straightforward, we generate testing data with predetermined +% The verification is straightforward: we generate test data with predetermined % flicker coefficients and verify the flicker assessment process by comparing -% the derived flicker coefficient with the predetermined ones. The testing data -% generated will be measured voltage and current time series (|u_m| and |i_m|), -% and the user can use |gen_test_data()| to generate test data for five different -% scenarios: (1) pure sine wave (coef=0) and scenarios (2)-(5) as stated Annex -% B.3.2 - B.3.5, and test them one by one. Here we illustrate the usage of gen_test_data -% by generating distorted u_m(t) with multiple zero crossings (see Annex B.3.3 -% for details). +% the derived flicker coefficients with the predetermined ones. The generated +% test data are measured voltage and current time series (|u_m| and |i_m|). The +% user can use |gen_test_data()| to generate test data for five scenarios: (1) a +% pure sine wave (coefficient = 0) and (2)-(5) the scenarios in Annex B.3.2 to +% B.3.5, and test them one by one. Here we illustrate |gen_test_data()| by +% generating a distorted u_m(t) with multiple zero crossings (Annex B.3.3). % Set up proper input parameters: % scenario 3 with duration of 10s opt = 3; T=10; -% rated parameters: +% rated parameters: voltage [V], current [A], sample rate [Hz], apparent power [VA] Un=12e3; In=144; fs=50e3; Sr = 3e6; % other parameters for generating distorted data fg=60; SCR=20; fm=25; fv=0.5; % fv only needed for scenario B.3.4 -% According to IEC TS 61400-21-1 Table B.2: +% According to IEC 61400-21-1:2019 Table B.2: DeltaI_I = [4.763 5.726 7.640 9.488]; [i_m,u_m]=gen_test_data(Un,In,fg,fs,fm,fv,DeltaI_I,opt,T) %% -% The generated i_m(ntime, nphi_k) has 4 time-series, corresponding to phi_k -% = 30, 50, 70, 85. While u_m is always the same, thus has only one time series. +% The generated i_m has four time series (columns), one for each impedance +% phase angle psi_k = 30, 50, 70, and 85 degrees. u_m is the same for all four +% angles, so it has a single time series. figure('Position', [100, 100, 1600, 600]); -hold off; plot(i_m.time, i_m.data(:,1)); xlim([0 1]); -xlabel("Time"); -ylabel("i_m"); +xlabel("Time [s]"); +ylabel("i_m [A]"); +title("Generated current, \psi_k = 30 degrees"); figure('Position', [100, 100, 1600, 600]); plot(u_m.time, u_m.data); xlim([0 0.1]); -xlabel("Time"); -ylabel("u_m"); +xlabel("Time [s]"); +ylabel("u_m [V]"); +title("Generated voltage with multiple zero crossings"); grid on; diff --git a/examples/power_example.mlx b/examples/power_example.mlx index 44f0e1517..e6d68eaee 100644 Binary files a/examples/power_example.mlx and b/examples/power_example.mlx differ diff --git a/mhkit/power/characteristics/ac_power_three_phase.m b/mhkit/power/characteristics/ac_power_three_phase.m index a002cd12d..9f1824e55 100644 --- a/mhkit/power/characteristics/ac_power_three_phase.m +++ b/mhkit/power/characteristics/ac_power_three_phase.m @@ -1,110 +1,114 @@ -function P=ac_power_three_phase(voltage,current,power_factor,varargin) +function P = ac_power_three_phase(voltage, current, power_factor, line_to_line) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Calculates the real power from three phase ac voltage and current. +% Calculate magnitude of active AC power from line to neutral voltage and current % % Parameters % ------------ -% voltage: Time series of all three measured voltages [V] -% Pandas data frame -% To make a pandas data frame from user supplied frequency and spectra -% use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(time,voltage) -% -% OR -% -% structure of form: -% voltage.voltage : matrix of all three phases -% -% voltage.time : time vector -% -% current: Time series of all three measured current [A] -% Pandas data frame -% To make a pandas data frame from user supplied frequency and spectra -% use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(time,current) -% -% OR -% -% structure of form: -% current.current : matrix of all three phases -% -% current.time : time vector -% -% power_factor : float -% power factor for the system -% -% line_to_line: bool (Optional) -% set true if the given voltage measurement is line_to_line +% voltage: structure or matrix +% voltage.voltage : Three-phase voltage measurements [V] (n_time x 3 matrix) +% Each row represents one time step, columns are phases A, B, C +% voltage.time : Time vector (n_time x 1) (if time series data) +% current: structure or matrix +% current.current : Three-phase current measurements [A] (n_time x 3 matrix) +% Each row represents one time step, columns are phases A, B, C +% current.time : Time vector (n_time x 1) (if time series data) +% power_factor: numeric scalar +% Power factor for the efficiency of the system [dimensionless] +% line_to_line: logical (optional) +% Set to true if voltage measurements are line-to-line (default: false) % % Returns % --------- -% P: Structure -% +% P: structure +% P.power : Magnitude of active AC power [W] +% P.time : Time vector % -% P.power [W] -% -% P.time -% -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -py.importlib.import_module('mhkit'); + arguments (Input) + voltage {mustBeA(voltage, ["struct", "numeric"])} + current {mustBeA(current, ["struct", "numeric"])} + power_factor (1,1) {mustBeNumeric} + line_to_line (1,1) {mustBeNumericOrLogical, mustBeMember(line_to_line, [0, 1])} = false + end -py.importlib.import_module('mhkit_python_utils'); + arguments (Output) + P struct + end -if (isa(voltage,'py.pandas.core.frame.DataFrame')~=1) - if (isstruct(voltage)==1) - x=size(voltage.voltage); - li=py.list(); - if x(2)==3 - for i = 1:x(2) - app=py.list(double(voltage.voltage(:,i))); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); + line_to_line = logical(line_to_line); - end - voltage=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(double(voltage.time(:,1)),li,int32(x(2))); - elseif x(2)==1 - ME = MException('MATLAB:ac_power_three_phase','Three lines of voltage measurements are required'); - throw(ME); - end + % Validate power factor is between 0 and 1 + if power_factor < 0 || power_factor > 1 + error('MHKiT:ac_power_three_phase:InvalidInput', 'power_factor must be between 0 and 1 (inclusive). Received: %.3f', power_factor); + end + % Extract data and time vectors + if isstruct(voltage) + % Validate input structures have required fields + if ~isfield(voltage, 'voltage') + error('MHKiT:ac_power_three_phase:InvalidInput', 'voltage structure must contain voltage field'); + end + if ~isfield(voltage, 'time') + error('MHKiT:ac_power_three_phase:InvalidInput', 'voltage structure must contain time field'); + end + voltage_data = voltage.voltage; + voltage_time = voltage.time; else - ME = MException('MATLAB:ac_power_three_phase','voltage needs to be a structure or Pandas dataframe, use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas to create one or a structure'); - throw(ME); + voltage_data = voltage; + voltage_time = (1:size(voltage, 1))'; % Default time vector end -end - -if (isa(current,'py.pandas.core.frame.DataFrame')~=1) - if (isstruct(current)==1) - x=size(current.current); - li=py.list(); - if x(2)==3 - for i = 1:x(2) - app=py.list(double(current.current(:,i))); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); - - end - current=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(double(current.time(:,1)),li,int32(x(2))); - elseif x(2)==1 - ME = MException('MATLAB:ac_power_three_phase','Three lines of voltage measurements are required'); - throw(ME); + + if isstruct(current) + % Validate input structures have required fields + if ~isfield(current, 'current') + error('MHKiT:ac_power_three_phase:InvalidInput', 'current structure must contain current field'); end - + if ~isfield(current, 'time') + error('MHKiT:ac_power_three_phase:InvalidInput', 'current structure must contain time field'); + end + current_data = current.current; + current_time = current.time; else - ME = MException('MATLAB:ac_power_three_phase','voltage needs to be a structure or Pandas dataframe, use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas to create one or a structure'); - throw(ME); + current_data = current; + current_time = (1:size(current, 1))'; % Default time vector end -end -if nargin == 4 - p_pd=py.mhkit.power.characteristics.ac_power_three_phase(voltage,current,power_factor,pyargs('line_to_line',varargin{1})); -elseif nargin ==3 - p_pd=py.mhkit.power.characteristics.ac_power_three_phase(voltage,current,power_factor); -else - ME = MException('MATLAB:ac_power_three_phase','incorrect number of input arguments'); - throw(ME); -end + % Validate voltage has three columns + if size(voltage_data, 2) ~= 3 + error('MHKiT:ac_power_three_phase:InvalidInput', 'voltage must have three columns for three-phase measurements'); + end + + % Validate current has three columns + if size(current_data, 2) ~= 3 + error('MHKiT:ac_power_three_phase:InvalidInput', 'current must have three columns for three-phase measurements'); + end + + % Validate dimensions match + if ~isequal(size(voltage_data), size(current_data)) + error('MHKiT:ac_power_three_phase:InvalidInput', 'voltage and current must have the same dimensions'); + end + + % Validate time vectors match + if ~isequal(voltage_time, current_time) + error('MHKiT:ac_power_three_phase:InvalidInput', 'voltage and current time vectors must match'); + end + + % After validation, use the time vector from voltage + time_vector = voltage_time; -P.power=double(py.array.array('d',py.numpy.nditer(p_pd.values))); -P.time=double(py.array.array('d',py.numpy.nditer(p_pd.index))); + % Instantaneous power summed across the three phases, as in MHKiT-Python + active_power = abs(sum(voltage_data .* current_data, 2)) * power_factor; + if line_to_line + active_power = active_power * sqrt(3); + end + + % Create output structure + P = struct(); + P.power = active_power; + P.time = time_vector; + +end diff --git a/mhkit/power/characteristics/dc_power.m b/mhkit/power/characteristics/dc_power.m index 353c731a1..ead1600e2 100644 --- a/mhkit/power/characteristics/dc_power.m +++ b/mhkit/power/characteristics/dc_power.m @@ -1,110 +1,80 @@ -function P=dc_power(voltage,current) +function power_dc = dc_power(voltage, current) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the real power from DC voltage and current. +% +% Calculate DC power from voltage and current % % Parameters % ------------ -% voltage: Time series of measured voltages [V] -% Pandas data frame -% To make a pandas data frame from user supplied frequency and spectra -% use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(time,voltage) -% -% OR -% -% structure of form: -% voltage.voltage : matrix or vector -% -% voltage.time : time vector -% -% current: Time series of current [A] -% Pandas data frame -% To make a pandas data frame from user supplied frequency and spectra -% use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(time,current) -% -% OR -% -% structure of form: -% current.current : matrix or vector -% -% current.time : time vector +% voltage: structure +% voltage.voltage : Matrix of voltage measurements [V] +% Each row is a time point, each column is a channel +% voltage.time : Time vector +% current: structure +% current.current : Matrix of current measurements [A] +% Each row is a time point, each column is a channel +% current.time : Time vector % % Returns % --------- -% P: Structure +% power_dc: structure +% power_dc.power : Matrix of calculated power for each channel [W] +% Same dimensions as input voltage/current +% power_dc.gross : Vector of gross power (sum of all channels) [W] +% power_dc.time : Time vector % -% -% P.power [W] -% -% P.gross: gross power from all lines [W] -% -% P.time -% -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -py.importlib.import_module('mhkit'); - -py.importlib.import_module('mhkit_python_utils'); - -if (isa(voltage,'py.pandas.core.frame.DataFrame')~=1) - if (isstruct(voltage)==1) - x=size(voltage.voltage); - li=py.list(); - if x(2)> 1 - for i = 1:x(2) - app=py.list(double(voltage.voltage(:,i))); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); - - end - voltage=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(double(voltage.time(:,1)),li,int32(x(2))); - elseif x(2)==1 - voltage=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(double(voltage.time(:,1)),voltage.voltage,int32(x(2))); - end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - else - ME = MException('MATLAB:dc_power','voltage needs to be a structure or Pandas dataframe, use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas to create one or create a structure'); - throw(ME); + arguments (Input) + voltage struct + current struct end -end -if (isa(current,'py.pandas.core.frame.DataFrame')~=1) - if (isstruct(current)==1) - x=size(current.current); - li=py.list(); - if x(2)==3 - for i = 1:x(2) - app=py.list(double(current.current(:,i))); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); - - end - current=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(double(current.time(:,1)),li,int32(x(2))); - elseif x(2)==1 - current=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(double(current.time(:,1)),current.current,int32(x(2))); - end - - else - ME = MException('MATLAB:dc_power','voltage needs to be a structure or Pandas dataframe, use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas to create one or create a structure'); - throw(ME); + arguments (Output) + power_dc struct end -end - -p_pd=py.mhkit.power.characteristics.dc_power(voltage,current); -vals=double(py.array.array('d',py.numpy.nditer(p_pd.values))); -sha=cell(p_pd.values.shape); -x=int64(sha{1,1}); -y=int64(sha{1,2}); -vals=reshape(vals,[x,y]); -si=size(vals); -% for i=1:si(2) -% if i == si(2) -% P.gross=vals(:,i); -% else -% P.power{i}=vals(:,i); -% %P.power{i}=[P.power{i}]; -% end -% end -P.power=vals(:,1:end-1); -P.gross=vals(:,end); -P.time=double(py.array.array('d',py.numpy.nditer(p_pd.index))).'; + + % Validate input structures have required fields + if ~isfield(voltage, 'voltage') + error('MHKiT:dc_power:InvalidInput', 'voltage structure must contain voltage field'); + end + if ~isfield(voltage, 'time') + error('MHKiT:dc_power:InvalidInput', 'voltage structure must contain time field'); + end + if ~isfield(current, 'current') + error('MHKiT:dc_power:InvalidInput', 'current structure must contain current field'); + end + if ~isfield(current, 'time') + error('MHKiT:dc_power:InvalidInput', 'current structure must contain time field'); + end + + % Extract data matrices + voltage_data = voltage.voltage; + current_data = current.current; + voltage_time = voltage.time; + current_time = current.time; + + % Validate dimensions match + if ~isequal(size(voltage_data), size(current_data)) + error('MHKiT:dc_power:InvalidInput', 'voltage and current must have the same dimensions'); + end + + % Validate time vectors match + if ~isequal(voltage_time, current_time) + error('MHKiT:dc_power:InvalidInput', 'Time vectors must match between voltage and current structures'); + end + + % Calculate power for each channel (element-wise multiplication) + power_matrix = voltage_data .* current_data; + + % Calculate gross power (sum across channels - dimension 2) + % Handle NaN values by skipping them in the sum + gross_power_vector = sum(power_matrix, 2, 'omitnan'); + + % Create output structure + power_dc = struct(); + power_dc.power = power_matrix; + power_dc.gross = gross_power_vector; + power_dc.time = voltage_time; +end diff --git a/mhkit/power/characteristics/instantaneous_frequency.m b/mhkit/power/characteristics/instantaneous_frequency.m index 2458607af..80c223150 100644 --- a/mhkit/power/characteristics/instantaneous_frequency.m +++ b/mhkit/power/characteristics/instantaneous_frequency.m @@ -1,51 +1,97 @@ -function frequency=instantaneous_frequency(voltage) +function frequency = instantaneous_frequency(voltage) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the instantaneous frequency of a measured voltage +% +% Calculate instantaneous frequency of measured voltage % % Parameters -% ----------- -% um: structure with handles- um.voltage and um.time -% measured voltage source (V) with each timeseries in its own column +% ------------ +% voltage : struct +% Measured voltage time series +% voltage.voltage : double [V] +% Measured voltage, n_samples x n_signals, one time series per column +% voltage.time : double [s] +% Time vector, n_samples x 1 % % Returns -% ------- -% frequency: structure -% frequency of the measured voltage with handles frequency.time and -% frequency.frequency -% +% --------- +% frequency : struct +% Instantaneous frequency of each voltage signal +% frequency.frequency : double [Hz] +% Instantaneous frequency, (n_samples - 1) x n_signals +% frequency.time : double [s] +% Time vector, (n_samples - 1) x 1. One element shorter than the +% input because the phase is differentiated. % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -py.importlib.import_module('mhkit'); -% py.importlib.import_module('numpy'); -py.importlib.import_module('mhkit_python_utils'); + arguments (Input) + voltage struct + end -time= voltage.time ; -data = voltage.voltage; -dname = 'voltage'; + arguments (Output) + frequency struct + end -dsize=size(data); + % Validate input structure has required fields + if ~isfield(voltage, 'voltage') + error('MHKiT:instantaneous_frequency:InvalidInput', 'voltage structure must contain voltage field'); + end + if ~isfield(voltage, 'time') + error('MHKiT:instantaneous_frequency:InvalidInput', 'voltage structure must contain time field'); + end -li=py.list(); -if dsize(2)>1 - for i = 1:dsize(2) - app=py.list(data(:,i)); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); + % Extract data from structure + voltage_data = voltage.voltage; + time_vector = voltage.time; - end - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(time(:,1),li,int32(dsize(2))); -elseif dsize(2)==1 - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(time,py.numpy.array(data),dsize(2)); -end + % Validate dimensions + if size(voltage_data, 1) ~= length(time_vector) + error('MHKiT:instantaneous_frequency:InvalidInput', 'voltage data rows must match time vector length'); + end + + % Get data dimensions + [num_samples, num_columns] = size(voltage_data); + + % Validate minimum data length for meaningful frequency calculation + if num_samples < 4 + error('MHKiT:instantaneous_frequency:InvalidInput', 'voltage data must have at least 4 samples for frequency calculation'); + end -harmonics_pd = py.mhkit.power.characteristics.instantaneous_frequency(data_pd); -vals=double(py.array.array('d',py.numpy.nditer(harmonics_pd.values))); -sha=cell(harmonics_pd.values.shape); -x=int64(sha{1,1}); -y=int64(sha{1,2}); -vals=reshape(vals,[x,y]); + % Warn the user if the sample interval varies by more than the tolerance. + % The phase derivative below divides by each local dt, but the FFT-based + % Hilbert transform that produces the phase assumes uniform sampling. + % Irregular intervals corrupt the phase itself before dt is applied, so + % using the local dt cannot correct for them. + sample_rate = mhkit_validate_sample_rate_hz(time_vector); + if ~sample_rate.is_uniform + warning('MHKiT:instantaneous_frequency:SampleRateVariation', ... + ['The sample interval of this signal varies by more than %g%% from the mean. ', ... + 'The FFT-based Hilbert transform assumes uniform sampling, so the ', ... + 'instantaneous phase, and therefore the instantaneous frequency, is ', ... + 'likely to be inaccurate. ', ... + 'Mean sample rate: %g Hz, max: %g Hz, min: %g Hz, standard deviation: %g Hz'], ... + sample_rate.tolerance_percent, sample_rate.mean_sample_rate_hz, ... + sample_rate.max_sample_rate_hz, sample_rate.min_sample_rate_hz, ... + sample_rate.std_sample_rate_hz); + end -frequency.frequency = vals; -frequency.time = time(2:end); + % Calculate time differences for frequency calculation + time_diff = diff(time_vector(:)); % Ensure column vector + + % Analytic signal of every column + analytic_signal = mhkit_hilbert(voltage_data); + + % Instantaneous phase with 2*pi discontinuities removed + unwrapped_phase = unwrap(angle(analytic_signal), [], 1); + + % Instantaneous frequency + frequency_data = diff(unwrapped_phase, 1, 1) ./ (2.0 * pi * time_diff); + + % Create output structure + frequency = struct(); + frequency.frequency = frequency_data; + frequency.time = time_vector(2:end); + +end diff --git a/mhkit/power/quality/calc_Rfic_Lfic.m b/mhkit/power/quality/calc_Rfic_Lfic.m index 44a6ee4ea..f309e3bd1 100644 --- a/mhkit/power/quality/calc_Rfic_Lfic.m +++ b/mhkit/power/quality/calc_Rfic_Lfic.m @@ -44,6 +44,18 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + arguments (Input) + Sr (1,1) {mustBeNumeric} + SCR (1,1) {mustBeNumeric} + Un (1,1) {mustBeNumeric} + fg (1,1) {mustBeNumeric} + end + + arguments (Output) + Rfic {mustBeNumeric} + Lfic {mustBeNumeric} + end + % impedance angles in degrees Phi_k = [30,50,70,85]; S_kfic = SCR*Sr; diff --git a/mhkit/power/quality/calc_electrical_angle.m b/mhkit/power/quality/calc_electrical_angle.m index ab6e5d70e..c267490eb 100644 --- a/mhkit/power/quality/calc_electrical_angle.m +++ b/mhkit/power/quality/calc_electrical_angle.m @@ -25,11 +25,18 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + arguments (Input) + freq struct + alpha0 (1,1) {mustBeNumeric} + end + + arguments (Output) + alpha_m {mustBeNumeric} + end + % check input: if ~isfield(freq,'time') || ~isfield(freq, 'data') - ME = MException('MATLAB:calc_electrail_angle',... - 'invalid handles in structure, must contain x.data & x.time'); - throw(ME); + error('MHKiT:calc_electrical_angle:InvalidInput', 'invalid handles in structure, must contain x.data & x.time'); end % IECTS-62600-30 Eq(3) diff --git a/mhkit/power/quality/calc_flicker_coefficient.m b/mhkit/power/quality/calc_flicker_coefficient.m index 3e71bfa55..f082b9ba4 100644 --- a/mhkit/power/quality/calc_flicker_coefficient.m +++ b/mhkit/power/quality/calc_flicker_coefficient.m @@ -30,6 +30,16 @@ % c(Phi_k) = P_stfic*S_kfic/Sr %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + arguments (Input) + P_stfic {mustBeNumeric} + S_kfic {mustBeNumeric} + Sr (1,1) {mustBeNumeric} + end + + arguments (Output) + coef_flicker {mustBeNumeric} + end + % IECTS62600-30 Eq (6): coef_flicker = P_stfic*S_kfic/Sr; end diff --git a/mhkit/power/quality/calc_fundamental_freq.m b/mhkit/power/quality/calc_fundamental_freq.m index ab8892dba..a54e88071 100644 --- a/mhkit/power/quality/calc_fundamental_freq.m +++ b/mhkit/power/quality/calc_fundamental_freq.m @@ -39,11 +39,20 @@ % Step3. Interpolate estimated freq back to u_m.time %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + arguments (Input) + u_m struct + method {mustBeTextScalar} + methodopts cell + end + + arguments (Output) + alpha0 {mustBeNumeric} + freq struct + end + % check input: if ~isfield(u_m,'time') || ~isfield(u_m, 'data') - ME = MException('MATLAB:calc_fundamental_freq',... - 'invalid handles in structure, must contain x.data & x.time'); - throw(ME); + error('MHKiT:calc_fundamental_freq:InvalidInput', 'invalid handles in structure, must contain x.data & x.time'); end % 1. make sure time starts at 0.0 u_m.time = u_m.time-u_m.time(1); diff --git a/mhkit/power/quality/calc_ideal_voltage.m b/mhkit/power/quality/calc_ideal_voltage.m index a94a7ecfb..738325dfb 100644 --- a/mhkit/power/quality/calc_ideal_voltage.m +++ b/mhkit/power/quality/calc_ideal_voltage.m @@ -24,6 +24,16 @@ % u0(t) = sqrt(2/3)*Un*sin(alpha_m(t)). % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + arguments (Input) + Un (1,1) {mustBeNumeric} + alpha_m {mustBeNumeric} + end + + arguments (Output) + u0 {mustBeNumeric} + end + %IECTS 62600-30(ed1.0) Eq (2) u0 = sqrt(2/3)*Un*sin(alpha_m); end diff --git a/mhkit/power/quality/calc_shortterm_flicker_severity.m b/mhkit/power/quality/calc_shortterm_flicker_severity.m index dd6ce47d9..f38ea6733 100644 --- a/mhkit/power/quality/calc_shortterm_flicker_severity.m +++ b/mhkit/power/quality/calc_shortterm_flicker_severity.m @@ -40,14 +40,20 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + arguments (Input) + P struct + end + + arguments (Output) + Pst {mustBeNumeric} + end + % check input: names_need = {'p30','p50','p80','p6','p8','p10','p13','p17',... 'p2p2','p3','p4','p0p7','p1','p1p5','p0p1'}; for i=1:length(names_need) if ~isfield(P,names_need{i}) - ME = MException('MATLAB:calc_shortterm_flicker_severity',... - 'invalid handles in structure, must contain x.data & x.time'); - throw(ME); + error('MHKiT:calc_shortterm_flicker_severity:InvalidInput', 'P structure must contain field %s', names_need{i}); end end P_50s = (P.p30+P.p50+P.p80)/3.; diff --git a/mhkit/power/quality/calc_simulated_voltage.m b/mhkit/power/quality/calc_simulated_voltage.m index 3b9112bd6..4e9d23b84 100644 --- a/mhkit/power/quality/calc_simulated_voltage.m +++ b/mhkit/power/quality/calc_simulated_voltage.m @@ -29,11 +29,20 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + arguments (Input) + u0 {mustBeNumeric} + i_m struct + Rfic {mustBeNumeric} + Lfic {mustBeNumeric} + end + + arguments (Output) + u_fic {mustBeNumeric} + end + % check input: if ~isfield(i_m,'time') || ~isfield(i_m, 'data') - ME = MException('MATLAB:calc_simulated_voltage',... - 'invalid handles in structure, must contain x.data & x.time'); - throw(ME); + error('MHKiT:calc_simulated_voltage:InvalidInput', 'invalid handles in structure, must contain x.data & x.time'); end % 1. calculate dim(t)/dt using diff() delta_im = diff(i_m.data);%reshape(,[],4); diff --git a/mhkit/power/quality/flicker_ufic_workflow.m b/mhkit/power/quality/flicker_ufic_workflow.m index 9f55e6d4b..933c06120 100644 --- a/mhkit/power/quality/flicker_ufic_workflow.m +++ b/mhkit/power/quality/flicker_ufic_workflow.m @@ -2,9 +2,9 @@ u_m,i_m,method,methodopts) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Conduct flicker assessment according to IECTS62600-30 and -% IECTS61400-21-1. -% A workflow of MATLAB functions up to calculation of u_fic. +% Conduct flicker assessment according to IEC TS 62600-30 and +% IEC 61400-21-1:2019. +% A workflow of MATLAB functions up to calculation of u_fic. % % Parameters % ----------- @@ -74,12 +74,25 @@ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +arguments (Input) + Sr (1,1) {mustBeNumeric} + Un (1,1) {mustBeNumeric} + SCR (1,1) {mustBeNumeric} + fg (1,1) {mustBeNumeric} + u_m struct + i_m struct + method {mustBeTextScalar} + methodopts cell +end + +arguments (Output) + out struct +end + % check input: if ~isfield(u_m,'time') || ~isfield(u_m, 'data') ||... ~isfield(i_m,'time') || ~isfield(i_m, 'data') - ME = MException('MATLAB:flicker_ufic',... - 'invalid handles in structure, must contain x.data & x.time'); - throw(ME); + error('MHKiT:flicker_ufic_workflow:InvalidInput', 'invalid handles in structure, must contain x.data & x.time'); end %% Step 1. Construct the fictitious grid: [Rfic,Lfic] = calc_Rfic_Lfic(Sr, SCR, Un, fg); diff --git a/mhkit/power/quality/gen_test_data.m b/mhkit/power/quality/gen_test_data.m index fbf2f6bad..c55336918 100644 --- a/mhkit/power/quality/gen_test_data.m +++ b/mhkit/power/quality/gen_test_data.m @@ -1,8 +1,8 @@ function [i_m,u_m]=gen_test_data(Un,In,fg,fs,fm,fv,DeltaI_I,opt,T) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Generate test data according to IEC61400-21-1 Annex B.3 to be used for -% the verification test of the measurement procedure for flicker. +% Generate test data according to IEC 61400-21-1:2019 Annex B.3 to be used for +% the verification test of the measurement procedure for flicker. % % Parameters % ----------- @@ -21,7 +21,7 @@ % into 60 different cases, with fν in [0.5:0.5:30]. % DeltaI_I: double array (4) % Relative current changes (%) according to Table B.2 and Table B.3 -% in IECTS61400-21-1. +% in IEC 61400-21-1:2019. % opt: int % Option number used to indicate different test datasets. % T: double @@ -49,6 +49,23 @@ % B.3.5 Slow frequency changes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + arguments (Input) + Un (1,1) {mustBeNumeric} + In (1,1) {mustBeNumeric} + fg (1,1) {mustBeNumeric} + fs (1,1) {mustBeNumeric} + fm (1,1) {mustBeNumeric} + fv (1,1) {mustBeNumeric} + DeltaI_I {mustBeNumeric} + opt (1,1) {mustBeNumeric} + T (1,1) {mustBeNumeric} + end + + arguments (Output) + i_m struct + u_m struct + end + %% time: t = 0:1/fs:T-1/fs; %600=10min t = t'; @@ -67,9 +84,7 @@ sin(2*pi*cumtrapz(t,func)); %plot(t,func);xlim([0 T]);grid on; else - ME = MException('MATLAB:gen_test_data',['option ' ... - 'does not exist']); - throw(ME); + error('MHKiT:gen_test_data:InvalidInput', 'option does not exist'); end %% u_m u_m = struct();u_m.time = t; @@ -96,9 +111,7 @@ sin(2*pi*cumtrapz(t,func)); % plot(t,u_m.data);xlim([0 1]);grid; else - ME = MException('MATLAB:gen_test_data',['option ' ... - 'does not exist']); - throw(ME); + error('MHKiT:gen_test_data:InvalidInput', 'option does not exist'); end end diff --git a/mhkit/power/quality/harmonic_subgroups.m b/mhkit/power/quality/harmonic_subgroups.m index 8229c644c..5718fef95 100644 --- a/mhkit/power/quality/harmonic_subgroups.m +++ b/mhkit/power/quality/harmonic_subgroups.m @@ -1,47 +1,95 @@ -function harmonic_subgroups=harmonic_subgroups(harmonics,grid_freq) +function harmonic_subgroups = harmonic_subgroups(harmonics, grid_freq) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the harmonic subgroups from the harmonics of current % -% Parameters -% ----------- -% harmonics: structure with handles- harmonics.amplitude and harmonics.harmonic -% harmonic amplitude with each timeseries in its own column +% Calculate the harmonic subgroups based on IEC 61000-4-7 % -% grid_freq: int -% value indicating if the power supply is 50 or 60 Hz. Valid input are 50 and 60 +% Parameters +% ------------ +% harmonics: structure +% harmonics.amplitude : Harmonic amplitude data with each timeseries in its own column +% harmonics.harmonic : Harmonic frequency vector [Hz] +% grid_freq: numeric +% Value indicating if the power supply is 50 or 60 Hz. Valid inputs are 50 and 60 % % Returns -% ------- -% harmonic_subgroups: structure with handles harmonic_subgroups.amplitude and -% harmonic_subgroups.harmonic +% --------- +% harmonic_subgroups: structure +% harmonic_subgroups.amplitude : Harmonic subgroup amplitudes +% harmonic_subgroups.harmonic : Target harmonic frequencies [Hz] % -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -data = harmonics.amplitude; + arguments (Input) + harmonics struct + grid_freq (1,1) {mustBeNumeric, mustBeMember(grid_freq, [50, 60])} + end -dsize=size(data); -li=py.list(); -if dsize(2)>1 - for i = 1:dsize(2) - app=py.list(data(:,i)); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); + arguments (Output) + harmonic_subgroups struct + end + + % Validate input structures have required fields + if ~isfield(harmonics, 'amplitude') + error('MHKiT:harmonic_subgroups:InvalidInput', 'harmonics structure must contain amplitude field'); + end + if ~isfield(harmonics, 'harmonic') + error('MHKiT:harmonic_subgroups:InvalidInput', 'harmonics structure must contain harmonic field'); + end + + % Extract data + harmonic_data = harmonics.amplitude; + harmonic_freq = harmonics.harmonic; + + % Validate dimensions + if size(harmonic_data, 1) ~= length(harmonic_freq) + error('MHKiT:harmonic_subgroups:InvalidInput', 'Number of harmonic frequencies must match the number of rows in amplitude data'); + end + + % Create target frequency array based on grid frequency + if grid_freq == 60 + target_hz = (0:60:3000)'; + else % grid_freq == 50 + target_hz = (0:50:2500)'; + end + + % Sort harmonics by frequency for proper indexing + [sorted_freq, sort_idx] = sort(harmonic_freq); + sorted_data = harmonic_data(sort_idx, :); + + % Get number of columns (signals) + num_signals = size(harmonic_data, 2); + num_targets = length(target_hz); + + % Initialize output array + subgroup_amplitudes = zeros(num_targets, num_signals); + + % Calculate harmonic subgroups + for i = 1:num_targets + current_target = target_hz(i); + + % Find nearest harmonic frequency index + [~, nearest_idx] = min(abs(sorted_freq - current_target)); + + % Define indices for the three adjacent harmonics (with bounds checking) + idx_prev = max(1, nearest_idx - 1); + idx_curr = nearest_idx; + idx_next = min(length(sorted_freq), nearest_idx + 1); + + % Calculate RMS sum for each signal column + for j = 1:num_signals + prev_val = sorted_data(idx_prev, j); + curr_val = sorted_data(idx_curr, j); + next_val = sorted_data(idx_next, j); + + % Calculate subgroup as RMS sum of three adjacent harmonics + subgroup_amplitudes(i, j) = sqrt(prev_val^2 + curr_val^2 + next_val^2); + end + end + + % Create output structure + harmonic_subgroups = struct(); + harmonic_subgroups.amplitude = subgroup_amplitudes; + harmonic_subgroups.harmonic = target_hz; - end - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(harmonics.harmonic,li,int32(dsize(2))); -elseif dsize(2)==1 - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(harmonics.harmonic,py.numpy.array(data),int32(dsize(2))); end - -sg_pd = py.mhkit.power.quality.harmonic_subgroups(data_pd,grid_freq); - -vals=double(py.array.array('d',py.numpy.nditer(sg_pd.values))); -sha=cell(sg_pd.values.shape); -x=int64(sha{1,1}); -y=int64(sha{1,2}); -vals=reshape(vals,[y,x]).'; - - -harmonic_subgroups.amplitude=vals; -harmonic_subgroups.harmonic = double(py.array.array('d',py.numpy.nditer(sg_pd.index))); diff --git a/mhkit/power/quality/harmonics.m b/mhkit/power/quality/harmonics.m index c3dd6feda..321a76034 100644 --- a/mhkit/power/quality/harmonics.m +++ b/mhkit/power/quality/harmonics.m @@ -1,70 +1,152 @@ -function harmonics=harmonics(x,freq,grid_freq) +function harmonics_result = harmonics(input_data, data_sample_rate_hz, grid_freq_hz, options) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the harmonics from time series of voltage or current based on IEC 61000-4-7. +% +% Calculate the harmonics from time series of voltage or current based on IEC 61000-4-7. % % Parameters -% ----------- -% x: structure with x.time and x.current or x.voltage as values -% Time-series of voltage [V] or current [A] +% ------------ +% input_data: structure +% input_data.current : Current time series [A], n_samples x n_signals +% (used if present) +% input_data.voltage : Voltage time series [V], n_samples x n_signals +% (used if current is not present) +% input_data.time : Time vector [s], n_samples x 1 (required) % -% freq: double -% Frequency of the time-series data [Hz] +% data_sample_rate_hz: double +% Sample rate of the time series [Hz]. IEC TS 62600-30 clause 7.1.4 +% requires at least 20 kHz per channel for harmonic measurements. % -% grid_freq: int -% Value indicating if the power supply is 50 or 60 Hz. Options = 50 or 60 +% grid_freq_hz: double +% Nominal grid frequency [Hz]. Options = 50 or 60 % -% Returns -% ------- -% harmonics: structure -% harmonic amplitude and frequency of the time-series data +% options.tolerance_percent: double (optional, default = 1) +% Allowed deviation [%] between the sample rate implied by +% input_data.time and data_sample_rate_hz. See mhkit_validate_sample_rate_hz. % +% Returns +% --------- +% harmonics_result: structure +% harmonics_result.amplitude : Amplitude A of the sinusoid A*sin(2*pi*f*t) +% at each frequency f, n_freqs x n_signals +% [A] for current, [V] for voltage. +% The 0 Hz row is abs(mean) of the signal. +% harmonics_result.harmonic : Frequencies [Hz], n_freqs x 1 +% = 0:5:((max_harmonic + 1) * grid_freq_hz - 5) +% harmonics_result.type : 'current' or 'voltage' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -py.importlib.import_module('mhkit'); -% py.importlib.import_module('numpy'); -py.importlib.import_module('mhkit_python_utils'); - -time= x.time ; - -if isfield(x, 'current') - data = x.current; - dname = 'current'; -elseif isfield(x,'voltage') - data = x.voltage; - dname = 'voltage'; -else - ME = MException('MATLAB:harmonics','invalid handles in structure, must contain x.current or x.voltage'); - throw(ME); -end + arguments (Input) + input_data struct + data_sample_rate_hz (1,1) double {mustBePositive} + grid_freq_hz (1,1) double {mustBeMember(grid_freq_hz, [50, 60])} + options.tolerance_percent (1,1) double {mustBePositive} = 1 + end -dsize=size(data); + arguments (Output) + harmonics_result struct + end -li=py.list(); -if dsize(2)>1 - for i = 1:dsize(2) - app=py.list(data(:,i)); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); + % Frequency spacing, IEC 61000-4-7 clause 3.4.1 NOTE 2 + freq_spacing_hz = 5; - end - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(time(:,1),li,int32(dsize(2))); -elseif dsize(2)==1 - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(time,py.numpy.array(data),dsize(2)); -end + % Highest harmonic order required by IEC TS 62600-30 clause 7.3 + max_harmonic = 50; + + % Validate input data structure + if ~isfield(input_data, 'time') + error('MHKiT:harmonics:InvalidInput', 'input_data structure must contain time field'); + end + + if isfield(input_data, 'current') + signal_data = input_data.current; + signal_type = 'current'; + elseif isfield(input_data, 'voltage') + signal_data = input_data.voltage; + signal_type = 'voltage'; + else + error('MHKiT:harmonics:InvalidInput', 'input_data structure must contain either current or voltage field'); + end + + if ~isnumeric(signal_data) || ~isreal(signal_data) + error('MHKiT:harmonics:InvalidInput', '%s data must be a real numeric array', signal_type); + end -harmonics_pd = py.mhkit.power.quality.harmonics(data_pd,freq,grid_freq); + % Validate time vector sample rate + sample_rate_validation = mhkit_validate_sample_rate_hz(input_data.time, data_sample_rate_hz, ... + 'tolerance_percent', options.tolerance_percent); -vals=double(py.array.array('d',py.numpy.nditer(harmonics_pd.values))); -sha=cell(harmonics_pd.values.shape); -x=int64(sha{1,1}); -y=int64(sha{1,2}); -vals=reshape(vals,[x,y]); + if ~sample_rate_validation.pass + error('MHKiT:harmonics:InvalidInput', ... + ['Time vector sample rate validation failed.\n' ... + 'Expected: %.2f Hz, Observed: %.2f Hz (median), Deviation: %.2f%%, Tolerance: %.2f%%\n' ... + 'Time format detected: %s'], ... + data_sample_rate_hz, sample_rate_validation.median_sample_rate_hz, ... + sample_rate_validation.deviation_percent, sample_rate_validation.tolerance_percent, ... + sample_rate_validation.time_format); + end + % Validate dimensions of every signal field against the time vector + num_time_samples = length(input_data.time); + field_names = fieldnames(input_data); + for field_idx = 1:length(field_names) + field_name = field_names{field_idx}; + if strcmp(field_name, 'time') + continue; + end + field_rows = size(input_data.(field_name), 1); + if field_rows ~= num_time_samples + error('MHKiT:harmonics:InvalidInput', '%s data rows (%d) must match time vector length (%d)', ... + field_name, field_rows, num_time_samples); + end + end -harmonics.amplitude=vals; -harmonics.harmonic = double(py.array.array('d',py.numpy.nditer(harmonics_pd.index))); -harmonics.type = dname; + num_signal_columns = size(signal_data, 2); + % Define frequencies [Hz] for harmonic amplitude reporting: 0, 5, 10, ... up to + % just below harmonic order max_harmonic + 1 (3055 Hz for a 60 Hz grid, + % 2545 Hz for a 50 Hz grid). + % + % 5 Hz spacing matches IEC 61000-4-7, which analyzes 200 ms windows + % (1 / 0.2 s = 5 Hz). + % + % Frequencies above harmonic order max_harmonic are included because + % harmonic_subgroups and interharmonics use them for order max_harmonic. + max_freq_hz = (max_harmonic + 1) * grid_freq_hz - freq_spacing_hz; + harmonic_freq_grid_hz = (0:freq_spacing_hz:max_freq_hz)'; + num_freqs = length(harmonic_freq_grid_hz); + % FFT amplitude of every column + fft_amplitude = abs(fft(signal_data, [], 1)); + + % Map each output frequency to the nearest FFT bin: bin k has frequency k * fs / n + exact_bin_index = harmonic_freq_grid_hz * num_time_samples / data_sample_rate_hz; + closest_bin_idx = round(exact_bin_index) + 1; % +1 for MATLAB 1-based indexing + nyquist_bin_idx = floor(num_time_samples / 2) + 1; + nyquist_freq_hz = data_sample_rate_hz / 2; + in_band = harmonic_freq_grid_hz <= nyquist_freq_hz & closest_bin_idx <= nyquist_bin_idx; + + if ~all(in_band) + warning('MHKiT:harmonics:AboveNyquist', ... + ['Frequencies from %.1f Hz upward are above the Nyquist frequency of the ' ... + '%.1f Hz sample rate and are set to zero. IEC TS 62600-30 clause 7.1.4 ' ... + 'requires at least 20 kHz per channel.'], ... + min(harmonic_freq_grid_hz(~in_band)), data_sample_rate_hz); + end + + harmonics_reindexed = zeros(num_freqs, num_signal_columns); + harmonics_reindexed(in_band, :) = fft_amplitude(closest_bin_idx(in_band), :); + + % Single-sided normalization: 2/n for all frequencies, 1/n for DC + harmonics_normalized = harmonics_reindexed / num_time_samples * 2; + harmonics_normalized(1, :) = harmonics_reindexed(1, :) / num_time_samples; + + % Create output structure + harmonics_result = struct(); + harmonics_result.amplitude = harmonics_normalized; + harmonics_result.harmonic = harmonic_freq_grid_hz; + harmonics_result.type = signal_type; + +end diff --git a/mhkit/power/quality/interharmonics.m b/mhkit/power/quality/interharmonics.m index 77f3c52e0..07802d02e 100644 --- a/mhkit/power/quality/interharmonics.m +++ b/mhkit/power/quality/interharmonics.m @@ -1,51 +1,116 @@ -function interharmonics=interharmonics(harmonics,grid_freq) +function interharmonics = interharmonics(harmonics, grid_freq) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the interharmonics from the harmonics of current based on IEC 61000-4-7. % -% Parameters -% ----------- -% harmonics: structure with handles- harmonics.amplitude and harmonics.harmonic -% Harmonic amplitude with each timeseries in its own column +% Calculate the interharmonics from the harmonic amplitudes of current % -% grid_freq: int -% Value indicating if the power supply is 50 or 60 Hz. Options = 50 or 60 +% Parameters +% ------------ +% harmonics: structure +% harmonics.amplitude : Harmonic amplitude data with each timeseries in its own column [dimensionless] +% harmonics.harmonic : Harmonic frequency vector corresponding to amplitude data [Hz] +% grid_freq: numeric scalar +% Value indicating if the power supply is 50 or 60 Hz. Options = 50 or 60 [Hz] % % Returns -% ------- -% interharmonics: structure with handles interharmonics.amplitude and -% interharmonics.harmonic -% +% --------- +% interharmonics: structure +% interharmonics.amplitude : Interharmonic group amplitudes [same units as input] +% interharmonics.harmonic : Frequency vector for interharmonic groups [Hz] % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -py.importlib.import_module('mhkit'); -% py.importlib.import_module('numpy'); -py.importlib.import_module('mhkit_python_utils'); + arguments (Input) + harmonics struct + grid_freq (1,1) {mustBeNumeric, mustBeMember(grid_freq, [50, 60])} + end -data = harmonics.amplitude; + arguments (Output) + interharmonics struct + end + + % Validate input structure has required fields + if ~isfield(harmonics, 'amplitude') + error('MHKiT:interharmonics:InvalidInput', 'harmonics structure must contain amplitude field'); + end + if ~isfield(harmonics, 'harmonic') + error('MHKiT:interharmonics:InvalidInput', 'harmonics structure must contain harmonic field'); + end + + % Extract data from structure + harmonics_amplitude = harmonics.amplitude; + harmonics_frequency = harmonics.harmonic; -dsize=size(data); -li=py.list(); -if dsize(2)>1 - for i = 1:dsize(2) - app=py.list(data(:,i)); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); + % Validate dimensions + if size(harmonics_amplitude, 1) ~= length(harmonics_frequency) + error('MHKiT:interharmonics:InvalidInput', 'harmonics.amplitude rows must match length of harmonics.harmonic'); + end + + % Create frequency grid based on grid frequency + if grid_freq == 60 + hz_grid = 0:60:3060; % MATLAB: 1-based, but this creates the correct range + hz_grid = hz_grid(1:end-1); % Remove last element, following MHKiT-Python convention + subset_size = 10; % Number of harmonics to include in each interharmonic group + elseif grid_freq == 50 + hz_grid = 0:50:2550; % MATLAB: 1-based, but this creates the correct range + hz_grid = hz_grid(1:end-1); % Remove last element, following MHKiT-Python convention + subset_size = 6; % Number of harmonics to include in each interharmonic group + end + + % Get number of columns in amplitude data + num_cols = size(harmonics_amplitude, 2); + num_frequencies = length(hz_grid); + + % Initialize output matrix + interharmonics_amplitude = ones(num_frequencies, num_cols); + + % Sort harmonics by frequency (equivalent to sort_index in Python) + [sorted_freq, sort_idx] = sort(harmonics_frequency); + sorted_amplitude = harmonics_amplitude(sort_idx, :); + + % Process each frequency in the grid + for i = 1:num_frequencies % MATLAB: 1-based indexing + current_freq = hz_grid(i); + + % Find nearest frequency index (equivalent to get_loc with method='nearest') + [~, nearest_idx] = min(abs(sorted_freq - current_freq)); + + % Process each column of amplitude data + for j = 1:num_cols % MATLAB: 1-based indexing + + % Calculate subset indices with bounds checking + % Python: harmonics[col].iloc[indn+1:indn+11] (60Hz) or [indn+1:indn+7] (50Hz) + % MATLAB: Need to add 1 for 1-based indexing and handle bounds + start_idx = nearest_idx + 1; % Python indn+1 becomes nearest_idx+1 in MATLAB + end_idx = start_idx + subset_size - 1; % Python slice end is exclusive, MATLAB inclusive + + % Ensure indices are within bounds + if start_idx > length(sorted_freq) + % If start is beyond array, use zeros + subset_squared = zeros(subset_size, 1); + elseif end_idx > length(sorted_freq) + % If end is beyond array, take what we can and pad with zeros + available_size = length(sorted_freq) - start_idx + 1; + if available_size > 0 + subset_data = sorted_amplitude(start_idx:end, j); + subset_squared = [subset_data.^2; zeros(subset_size - available_size, 1)]; + else + subset_squared = zeros(subset_size, 1); + end + else + % Normal case: extract subset and square + subset_data = sorted_amplitude(start_idx:end_idx, j); + subset_squared = subset_data.^2; % MATLAB: .^ for element-wise power + end + + % Calculate interharmonic as square root of sum of squares + interharmonics_amplitude(i, j) = sqrt(sum(subset_squared)); + end + end + + % Create output structure + interharmonics = struct(); + interharmonics.amplitude = interharmonics_amplitude; + interharmonics.harmonic = hz_grid'; % Convert to column vector to match input format - end - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(harmonics.harmonic,li,int32(dsize(2))); -elseif dsize(2)==1 - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(harmonics.harmonic,py.numpy.array(data),int32(dsize(2))); end - -interharmonics_pd = py.mhkit.power.quality.interharmonics(data_pd,grid_freq); -vals=double(py.array.array('d',py.numpy.nditer(interharmonics_pd.values))); -sha=cell(interharmonics_pd.values.shape); -x=int64(sha{1,1}); -y=int64(sha{1,2}); -vals=reshape(vals,[x,y]); - - -interharmonics.amplitude=vals; -interharmonics.harmonic = double(py.array.array('d',py.numpy.nditer(interharmonics_pd.index))); - diff --git a/mhkit/power/quality/total_harmonic_current_distortion.m b/mhkit/power/quality/total_harmonic_current_distortion.m index 518646ce5..8e63eb467 100644 --- a/mhkit/power/quality/total_harmonic_current_distortion.m +++ b/mhkit/power/quality/total_harmonic_current_distortion.m @@ -1,38 +1,95 @@ -function THCD=total_harmonic_current_distortion(harmonic_subgroups,rated_current) +function THCD = total_harmonic_current_distortion(harmonic_subgroups, rated_current) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the total harmonic current distortion (THC) based on IEC/TS 62600-30 % -% Parameters -% ----------- -% harmonic_subgroups: structure with handles- harmonic_subgroups.amplitude and harmonic_subgroups.harmonic -% Subgrouped current harmonics amplitude indexed by harmonic order +% Calculate the total harmonic current distortion (THC) based on IEC TS 62600-30 % -% rated_current: double -% Rated current of the energy device in Amps +% Parameters +% ------------ +% harmonic_subgroups: structure +% harmonic_subgroups.amplitude : Subgrouped current harmonics amplitude indexed by harmonic order [Amps] +% harmonic_subgroups.harmonic : Harmonic frequency order vector +% rated_current: double +% Rated current of the energy device [Amps] % % Returns -% ------- -% THCD: double -% Total harmonic current distortion -% +% --------- +% THCD: double +% Total harmonic current distortion [%] % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -data = harmonic_subgroups.amplitude; + arguments (Input) + harmonic_subgroups struct + rated_current (1,1) {mustBeNumeric, mustBePositive} + end -dsize=size(data); -li=py.list(); -if dsize(2)>1 - for i = 1:dsize(2) - app=py.list(data(:,i)); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); + arguments (Output) + THCD + end + + % Validate input structures have required fields + if ~isfield(harmonic_subgroups, 'amplitude') + error('MHKiT:total_harmonic_current_distortion:InvalidInput', 'harmonic_subgroups structure must contain amplitude field'); + end + if ~isfield(harmonic_subgroups, 'harmonic') + error('MHKiT:total_harmonic_current_distortion:InvalidInput', 'harmonic_subgroups structure must contain harmonic field'); + end + + % Extract amplitude data + harmonics_data = harmonic_subgroups.amplitude; + + % Validate amplitude data dimensions + if ~isnumeric(harmonics_data) + error('MHKiT:total_harmonic_current_distortion:InvalidInput', 'harmonic_subgroups.amplitude must be numeric'); + end + + % Check if we have enough harmonic data (need at least index 2 for fundamental and some harmonics) + if length(harmonics_data) < 3 + error('MHKiT:total_harmonic_current_distortion:InvalidInput', 'harmonic_subgroups.amplitude must contain at least 3 elements'); + end + + % Convert Python indexing to MATLAB indexing + % Python: harmonics_subgroup.iloc[2:50] means indices 2 through 49 (0-based) + % MATLAB: equivalent is indices 3 through 50 (1-based) + harmonic_start_idx = 3; % Python index 2 + 1 + harmonic_end_idx = min(50, length(harmonics_data)); % Python index 49 + 1, but limit to data length + + % Extract harmonic components (excluding fundamental) + % Python: harmonics_subgroup.iloc[2:50]**2 + % MATLAB: harmonics_data(3:end_idx).^2 (element-wise power) + if harmonic_end_idx >= harmonic_start_idx + harmonics_subset = harmonics_data(harmonic_start_idx:harmonic_end_idx); + else + error('MHKiT:total_harmonic_current_distortion:InvalidInput', 'Insufficient harmonic data for calculation'); + end + + % Square the harmonic amplitudes (element-wise operation) + harmonics_sq = harmonics_subset .^ 2; + + % Sum the squared harmonics + harmonics_sum = sum(harmonics_sq); + + % Get fundamental component + % Python: harmonics_subgroup.iloc[1] means index 1 (0-based) + % MATLAB: harmonics_data(2) means index 2 (1-based) + if length(harmonics_data) < 2 + error('MHKiT:total_harmonic_current_distortion:InvalidInput', 'harmonic_subgroups.amplitude must contain fundamental component at index 2'); + end + + fundamental_current = harmonics_data(2); % Python index 1 + 1 + + % Validate fundamental current is not zero + if fundamental_current == 0 + error('MHKiT:total_harmonic_current_distortion:InvalidInput', 'Fundamental current component cannot be zero'); + end + + % Calculate Total Harmonic Current Distortion + % Python: (np.sqrt(harmonics_sum)/harmonics_subgroup.iloc[1])*100 + % MATLAB: (sqrt(harmonics_sum) ./ fundamental_current) .* 100 + THCD = (sqrt(harmonics_sum) ./ fundamental_current) .* 100; + + % Ensure output is a scalar double + THCD = double(THCD); - end - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(harmonic_subgroups.harmonic,li,int32(dsize(2))); -elseif dsize(2)==1 - data_pd=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(harmonic_subgroups.harmonic,py.numpy.array(data),int32(dsize(2))); end - -thcd_pd = py.mhkit.power.quality.total_harmonic_current_distortion(data_pd); -THCD = double(py.array.array('d',py.numpy.nditer(thcd_pd.values))); diff --git a/mhkit/signal_processing/mhkit_hilbert.m b/mhkit/signal_processing/mhkit_hilbert.m new file mode 100644 index 000000000..b739d560c --- /dev/null +++ b/mhkit/signal_processing/mhkit_hilbert.m @@ -0,0 +1,58 @@ +function analytic_signal = mhkit_hilbert(x) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Compute the analytic signal of a real signal using the Hilbert transform +% +% Native MATLAB implementation that needs neither the Signal Processing +% Toolbox nor Python. It returns the analytic signal z = x + i*H[x], where +% H is the Hilbert transform, which is the same quantity returned by the +% Signal Processing Toolbox hilbert function and by scipy.signal.hilbert. +% The algorithm matches scipy.signal.hilbert in SciPy v1.16.2: +% https://github.com/scipy/scipy/blob/b1296b9b4393e251511fe8fdd3e58c22a1124899/scipy/signal/_signaltools.py#L2476-L2594 +% +% Parameters +% ------------ +% x : double +% Real, non-empty signal. A vector is transformed along its length. A +% matrix is transformed column by column, one time series per column. +% +% Returns +% --------- +% analytic_signal : complex double +% Analytic signal, same size as x. real(analytic_signal) equals x and +% imag(analytic_signal) is the Hilbert transform of x. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + arguments + x {mustBeNumeric, mustBeReal} + end + + % Transform a row vector along its length, like hilbert() and scipy + [x, was_row] = mhkit_standardize_user_input_to_column_vectors(double(x), ... + 'function_name', 'mhkit_hilbert'); + + n = size(x, 1); + spectrum = fft(x, [], 1); + + h = zeros(n, 1); + if mod(n, 2) == 0 + h([1, n / 2 + 1]) = 1; + h(2:n / 2) = 2; + else + h(1) = 1; + h(2:(n + 1) / 2) = 2; + end + + analytic_signal = ifft(spectrum .* h, [], 1); + + % Always return complex, as scipy does, even when the imaginary part is + % exactly zero (for example n = 1 or n = 2) + if isreal(analytic_signal) + analytic_signal = complex(analytic_signal); + end + + analytic_signal = mhkit_restore_column_vectors_to_user_input(analytic_signal, was_row); + +end diff --git a/mhkit/tests/Power_TestDevice.m b/mhkit/tests/Power_TestDevice.m index 850da458c..8b16256e9 100644 --- a/mhkit/tests/Power_TestDevice.m +++ b/mhkit/tests/Power_TestDevice.m @@ -66,6 +66,104 @@ function test_interharmonics_sine_wave(testCase) assertEqual(testCase,inter_harmonics.amplitude', interharmonic,'AbsTol',0.1) end + function test_harmonics_multi_harmonic_60hz(testCase) + % DC + fundamental + orders 3 and 5 + one interharmonic at 1.5 x + % the fundamental, sampled at 10 kHz for 1 s + fs = 10000; + frequency = 60; + samples = (0:fs-1)' / fs; + current = 0.3 ... + + 1.0 * sin(2*pi*frequency*samples) ... + + 0.10 * sin(2*pi*3*frequency*samples) ... + + 0.05 * sin(2*pi*5*frequency*samples) ... + + 0.02 * sin(2*pi*1.5*frequency*samples); + x = struct('current', current, 'time', samples); + + h = harmonics(x, fs, frequency); + + testCase.verifyEqual(size(h.amplitude), [612 1]); + testCase.verifyEqual(size(h.harmonic), [612 1]); + testCase.verifyEqual(h.type, 'current'); + + expected_hz = [0, 60, 180, 300, 90]; + expected_amp = [0.3, 1.0, 0.10, 0.05, 0.02]; + for k = 1:numel(expected_hz) + testCase.verifyEqual(h.amplitude(h.harmonic == expected_hz(k)), ... + expected_amp(k), 'AbsTol', 1e-3, ... + sprintf('Amplitude at %d Hz', expected_hz(k))); + end + other = ~ismember(h.harmonic, expected_hz); + testCase.verifyLessThan(max(abs(h.amplitude(other))), 1e-3); + end + + function test_harmonics_warns_above_nyquist(testCase) + % At fs = 1000 Hz the lines above 500 Hz cannot be resolved + fs = 1000; + t = (0:fs-1)' / fs; + x = struct('current', sin(2*pi*60*t), 'time', t); + + h = testCase.verifyWarning(@() harmonics(x, fs, 60), ... + 'MHKiT:harmonics:AboveNyquist'); + testCase.verifyEqual(numel(h.harmonic), 612); + testCase.verifyEqual(h.amplitude(h.harmonic > fs / 2), zeros(511, 1)); + testCase.verifyEqual(h.amplitude(h.harmonic == 60), 1, 'AbsTol', 1e-9); + end + + function test_harmonics_errors(testCase) + fs = 10000; + samples = (0:fs-1)' / fs; + current = sin(2*pi*60*samples); + + % (1) Missing time field + x_no_time = struct('current', current); + testCase.verifyError(@() harmonics(x_no_time, fs, 60), ?MException); + + % (2) Neither current nor voltage + x_no_signal = struct('time', samples); + testCase.verifyError(@() harmonics(x_no_signal, fs, 60), ?MException); + + % (3) Time vector at 9000 Hz passed with fs = 10000 Hz + t9000 = (0:8999)' / 9000; + x9000 = struct('current', sin(2*pi*60*t9000), 'time', t9000); + testCase.verifyError(@() harmonics(x9000, 10000, 60), ?MException); + end + + function test_mhkit_validate_sample_rate_hz(testCase) + % Pass case + time_s = (0:1/1000:1)'; + result = mhkit_validate_sample_rate_hz(time_s, 1000); + testCase.verifyTrue(result.pass); + testCase.verifyEqual(result.median_sample_rate_hz, 1000, 'RelTol', 1e-6); + testCase.verifyEqual(result.mean_sample_rate_hz, 1000, 'RelTol', 1e-6); + testCase.verifyTrue(result.is_uniform); + + % Expected rate omitted: defaults to the nominal rate of the record + result_nominal = mhkit_validate_sample_rate_hz(time_s); + testCase.verifyTrue(result_nominal.pass); + testCase.verifyTrue(result_nominal.is_uniform); + + % Outlier tolerance: one irregular interval must not fail the record + time_outlier = (0:1/1000:1)'; + time_outlier(500) = time_outlier(500) + 0.4 / 1000; + result = mhkit_validate_sample_rate_hz(time_outlier, 1000); + testCase.verifyTrue(result.pass); + testCase.verifyEqual(result.median_sample_rate_hz, 1000, 'RelTol', 1e-6); + testCase.verifyGreaterThan(result.max_deviation_percent, 1); + testCase.verifyEqual(result.time_format, "seconds elapsed"); + testCase.verifyFalse(result.is_uniform); + testCase.verifyEqual(result.max_interval_deviation_percent, 40, 'RelTol', 1e-2); + + % Fail case: 1000 Hz data against an expected 1100 Hz + result_fail = mhkit_validate_sample_rate_hz(time_s, 1100); + testCase.verifyFalse(result_fail.pass); + + % POSIX timestamp case + time_posix = 1.7e9 + (0:0.001:1)'; + result_posix = mhkit_validate_sample_rate_hz(time_posix, 1000); + testCase.verifyTrue(result_posix.pass); + testCase.verifyEqual(result_posix.time_format, "POSIX timestamp"); + end + function test_instfreq(testCase) frequency = 60; t = 600; @@ -124,13 +222,13 @@ function test_gen_test_data(testCase) i_m0 = readmatrix("../../examples/data/power/testdata/sinewave-pi6_im.txt"); testCase.verifyTrue(max(abs( ... (i_m.data(1:100,1)-i_m0(:,1))./i_m0(:,1) ... - ))<1e-10,opt); + ))<1e-10,string(opt)); testCase.verifyTrue(max(abs( ... (u_m.data(1:100)-u_m0)./u_m0)) ... - <1e-10,opt); + <1e-10,string(opt)); - % 2. opt=randi([2,5]) tests generated according to IECTS - opt = randi([2,5],1); + % 2. opt=2:5 tests generated according to IECTS + for opt = 2:5 if opt==2 SCR = 50; fm = 33.3; % TableB.3,fg=60, SCR=50 @@ -148,16 +246,17 @@ function test_gen_test_data(testCase) sprintf('../../examples/data/power/testdata/B.3.%i_um.txt',opt)); testCase.verifyTrue(max(abs( ... (i_m.data(1:100,1)-i_m0(:,1))./i_m0(:,1) ... - ))<1e-10,opt); + ))<1e-10,string(opt)); testCase.verifyTrue(max(abs((u_m.data(1:100)-u_m0)./u_m0))<1e-10,... - opt); + string(opt)); + end end function test_flicker_ufic_workflow(testCase) % u_m, i_m, Sr, Un, In, SCR, fg, & fs Sr = 3e6; Un=12e3; In=144; fg=60; fs=50e3;fv=0.5; %1. opt = 0, pure sine waves: - opt = 0; idx = randi([1,4],1); + opt = 0; fm=20; SCR=20; DeltaI_I = [0 0 0 0];% pure sine wave [i_m,u_m]=gen_test_data(Un,In,fg,fs,fm,fv,DeltaI_I,opt,10); method = 'ZCD'; methodopts = {}; @@ -173,11 +272,14 @@ function test_flicker_ufic_workflow(testCase) (out.alpha_m(1:100)-alpha_m0)./alpha_m0))<1e-10,string(opt)); testCase.verifyTrue(max(abs( ... (out.u0(1:100)-u00)./u00))<1e-10,string(opt)); - testCase.verifyTrue(max(abs(( ... - out.u_fic(1:100,idx)-u_fic0(:,idx))./u_fic0(:,idx)))<1e-10, ... - string(opt)); - % 2. opt=randi([2,5]) tests generated according to IECTS - opt = randi([2,5],1); idx = randi([1,4],1); + for idx = 1:4 + testCase.verifyTrue(max(abs(( ... + out.u_fic(1:100,idx)-u_fic0(:,idx))./u_fic0(:,idx)))<1e-10, ... + string(opt)); + end + % 2. opt=2,4,5 tests generated according to IECTS (opt=3 needs + % STFT, see test_flicker_ufic_workflow_stft) + for opt = [2, 4, 5] if opt==2 SCR = 50; fm = 33.3; % TableB.3,fg=60,SCR=50 @@ -189,14 +291,47 @@ function test_flicker_ufic_workflow(testCase) %DeltaI_I = [3.212 3.958 5.644 7.711];% fm=20 end [i_m,u_m]=gen_test_data(Un,In,fg,fs,fm,fv,DeltaI_I,opt,10); - if opt==3 - % B.3.3 Distorted um with multiple zero crossings - method = 'stft';methodopts = {... - 'Window',rectwin(int32(50000)),'OverlapLength',25000,... - 'FFTLength',50e3,'FrequencyRange','onesided'}; - else - method = 'ZCD'; methodopts = {}; + method = 'ZCD'; methodopts = {}; + out = flicker_ufic_workflow(... + Sr,Un,SCR,fg,u_m,i_m,method,methodopts); + freq0 = readmatrix( ... + sprintf('../../examples/data/power/testdata/B.3.%i_freq.txt',opt)); + alpha_m0 = readmatrix( ... + sprintf('../../examples/data/power/testdata/B.3.%i_alpham.txt',opt)); + u00 = readmatrix( ... + sprintf('../../examples/data/power/testdata/B.3.%i_u0.txt',opt)); + u_fic0 = readmatrix(... + sprintf('../../examples/data/power/testdata/B.3.%i_ufic.txt',opt)); + testCase.verifyTrue(max(abs( ... + (out.freq.data(1:100)-freq0)./freq0))<1e-10, string(opt)); + testCase.verifyTrue(max(abs( ... + (out.alpha_m(1:100)-alpha_m0)./alpha_m0))<1e-10, string(opt)); + testCase.verifyTrue(max(abs( ... + (out.u0(1:100)-u00)./u00))<1e-10, string(opt)); + for idx = 1:4 + testCase.verifyTrue(max(abs(( ... + out.u_fic(1:100,idx)-u_fic0(:,idx))./u_fic0(:,idx)))<1e-10, ... + string(opt)); + end end + end + + function test_flicker_ufic_workflow_stft(testCase) + % B.3.3 Distorted um with multiple zero crossings requires the + % STFT method, which needs the Signal Processing Toolbox + testCase.assumeTrue(exist('stft', 'file') ~= 0 && exist('rectwin', 'file') ~= 0, ... + 'Requires the Signal Processing Toolbox (stft, rectwin)'); + + % u_m, i_m, Sr, Un, In, SCR, fg, & fs + Sr = 3e6; Un=12e3; In=144; fg=60; fs=50e3;fv=0.5; + opt = 3; + SCR = 20; fm = 25; + % TableB.2,fg=60,SCR=20 + DeltaI_I = [4.763 5.726 7.640 9.488];% fm=25 + [i_m,u_m]=gen_test_data(Un,In,fg,fs,fm,fv,DeltaI_I,opt,10); + method = 'stft';methodopts = {... + 'Window',rectwin(int32(50000)),'OverlapLength',25000,... + 'FFTLength',50e3,'FrequencyRange','onesided'}; out = flicker_ufic_workflow(... Sr,Un,SCR,fg,u_m,i_m,method,methodopts); freq0 = readmatrix( ... @@ -213,9 +348,11 @@ function test_flicker_ufic_workflow(testCase) (out.alpha_m(1:100)-alpha_m0)./alpha_m0))<1e-10, string(opt)); testCase.verifyTrue(max(abs( ... (out.u0(1:100)-u00)./u00))<1e-10, string(opt)); - testCase.verifyTrue(max(abs(( ... - out.u_fic(1:100,idx)-u_fic0(:,idx))./u_fic0(:,idx)))<1e-10, ... - string(opt)); + for idx = 1:4 + testCase.verifyTrue(max(abs(( ... + out.u_fic(1:100,idx)-u_fic0(:,idx))./u_fic0(:,idx)))<1e-10, ... + string(opt)); + end end function test_calc_Rfic_Lfic(testCase) diff --git a/mhkit/tests/SignalProcessing_TestHilbert.m b/mhkit/tests/SignalProcessing_TestHilbert.m new file mode 100644 index 000000000..9775e39cb --- /dev/null +++ b/mhkit/tests/SignalProcessing_TestHilbert.m @@ -0,0 +1,70 @@ +classdef SignalProcessing_TestHilbert < matlab.unittest.TestCase + % Tests for mhkit_hilbert using exact analytic-signal identities. + % + % For a signal with an integer number of periods in the record: + % cos(w t) -> exp(i w t) (H[cos] = sin) + % sin(w t) -> -i exp(i w t) (H[sin] = -cos) + % and the Hilbert transform of a constant (DC) is zero. + + properties (TestParameter) + n_samples = {63, 64, 1000, 1001}; + end + + methods (Test) + function test_cosine_and_sine(testCase, n_samples) + k = 5; % whole periods in the record + t = (0:n_samples-1)' / n_samples; + phase = 2 * pi * k * t; + + testCase.verifyEqual(mhkit_hilbert(cos(phase)), exp(1i * phase), 'AbsTol', 1e-12); + testCase.verifyEqual(mhkit_hilbert(sin(phase)), -1i * exp(1i * phase), 'AbsTol', 1e-12); + end + + function test_dc_has_zero_hilbert_transform(testCase, n_samples) + x = 3.5 * ones(n_samples, 1); + z = mhkit_hilbert(x); + testCase.verifyEqual(real(z), x, 'AbsTol', 1e-12); + testCase.verifyEqual(imag(z), zeros(n_samples, 1), 'AbsTol', 1e-12); + end + + function test_real_part_is_input(testCase, n_samples) + x = randn(RandStream('twister', 'Seed', n_samples), n_samples, 1); + testCase.verifyEqual(real(mhkit_hilbert(x)), x, 'AbsTol', 1e-12); + end + + function test_matrix_is_transformed_column_wise(testCase) + n = 256; + t = (0:n-1)' / n; + x = [cos(2*pi*3*t), 2 * cos(2*pi*7*t), sin(2*pi*11*t)]; + z = mhkit_hilbert(x); + testCase.verifySize(z, [n, 3]); + for column = 1:3 + testCase.verifyEqual(z(:, column), mhkit_hilbert(x(:, column)), 'AbsTol', 1e-12); + end + end + + function test_row_vector_is_transformed_along_its_length(testCase) + n = 128; + phase = 2 * pi * 4 * (0:n-1) / n; + z = mhkit_hilbert(cos(phase)); + testCase.verifySize(z, [1, n]); + testCase.verifyEqual(z, exp(1i * phase), 'AbsTol', 1e-12); + end + + function test_integer_input(testCase) + x = int16([3; -7; 2; 9; -1; 0; 5]); + testCase.verifyEqual(mhkit_hilbert(x), mhkit_hilbert(double(x))); + end + + function test_short_inputs(testCase) + % n = 1 and n = 2 have no negative frequencies to remove + testCase.verifyEqual(mhkit_hilbert(4), complex(4, 0)); + testCase.verifyEqual(mhkit_hilbert([1; -2]), complex([1; -2], 0), 'AbsTol', 1e-15); + end + + function test_invalid_input_errors(testCase) + testCase.verifyError(@() mhkit_hilbert([1; 1i]), 'MATLAB:validators:mustBeReal'); + testCase.verifyError(@() mhkit_hilbert([]), 'MHKiT:mhkit_hilbert:InvalidInput'); + end + end +end diff --git a/mhkit/tests/examples/Examples_TestMatlabOnly.m b/mhkit/tests/examples/Examples_TestMatlabOnly.m index 8a9a48e18..59cf85ac5 100644 --- a/mhkit/tests/examples/Examples_TestMatlabOnly.m +++ b/mhkit/tests/examples/Examples_TestMatlabOnly.m @@ -16,6 +16,7 @@ 'adv_example', ... 'cdip_example', ... 'mooring_example', ... + 'power_example', ... 'strain_measurement_example', ... 'upcrossing_example', ... }; diff --git a/mhkit/tests/runMatlabOnlyTests.m b/mhkit/tests/runMatlabOnlyTests.m index d20cfd531..cc2a9a8dd 100644 --- a/mhkit/tests/runMatlabOnlyTests.m +++ b/mhkit/tests/runMatlabOnlyTests.m @@ -23,6 +23,8 @@ 'Dolfyn_Test_Rotate.m', ... 'Dolfyn_Test_VAP.m', ... 'Mooring_TestMooring.m', ... + 'Power_TestDevice.m', ... + 'SignalProcessing_TestHilbert.m', ... 'Tidal_TestIO.m', ... 'upcrossing_Test.m', ... }; diff --git a/mhkit/utils/mhkit_restore_column_vectors_to_user_input.m b/mhkit/utils/mhkit_restore_column_vectors_to_user_input.m new file mode 100644 index 000000000..9b1e3836c --- /dev/null +++ b/mhkit/utils/mhkit_restore_column_vectors_to_user_input.m @@ -0,0 +1,43 @@ +function data = mhkit_restore_column_vectors_to_user_input(data, was_row) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Restores column-oriented data to the caller's original vector orientation +% +% Sister function to mhkit_standardize_user_input_to_column_vectors. Only +% use this for shape-preserving functions, where the output has the same +% number of elements as the standardized input (e.g. surface_elevation). +% Spectral statistic functions (fewer output elements than input, e.g. +% frequency_moment, peak_period) always emit column-oriented output per +% MHKiT's vector convention and must not call this function. +% +% Parameters +% ------------ +% data : numeric vector or matrix +% Column-oriented output to restore +% was_row : logical +% Orientation flag returned by +% mhkit_standardize_user_input_to_column_vectors +% +% Returns +% --------- +% data : numeric vector or matrix +% Transposed back to a row vector if was_row is true, unchanged +% otherwise +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + data {mustBeNumeric} + was_row (1,1) logical +end + +arguments (Output) + data +end + +if was_row + data = data.'; +end + +end diff --git a/mhkit/utils/mhkit_standardize_user_input_to_column_vectors.m b/mhkit/utils/mhkit_standardize_user_input_to_column_vectors.m new file mode 100644 index 000000000..1f46ac087 --- /dev/null +++ b/mhkit/utils/mhkit_standardize_user_input_to_column_vectors.m @@ -0,0 +1,66 @@ +function [data, was_row] = mhkit_standardize_user_input_to_column_vectors(data, options) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Standardizes numeric vector input to MHKiT's column-vector convention +% +% MHKiT-MATLAB functions operate on column vectors (Nx1) internally, and +% on matrices (NxK) representing K column-oriented vectors (e.g. K wave +% spectra of length N, one per column). A caller may pass a row vector +% (1xN) instead; this function transposes it to a column and reports the +% original orientation so the caller's output can later be restored via +% mhkit_restore_column_vectors_to_user_input. Matrices (MxN, M>1 and +% N>1) are already column-oriented by convention and pass through +% unchanged. +% +% Parameters +% ------------ +% data : numeric vector or matrix +% User-supplied numeric input to standardize +% function_name : string (optional) +% Name of the calling function, used to scope any error this function +% raises (e.g. "MHKiT: