%% Estimation of electrochemical cell impedance
%  using CWT with complex Morlet mother wavelet
%
%  This .m file is intended to show how the impedance can be estimated
%  based on the procedure described by Debenjak et. el (2014). The paper
%  talks about PEM fuel cell impedance, however, the exact procedure can be
%  used on any kind of electrochemical cell.

%  In the .mat files are voltage and current measurements taken on a
%  laboratory-grade Li-S cell. Voltage and Current are measured in mV and
%  mA, respectively. Each .mat file contains measurements at specific PRBS
%  sampling frequency covering specific frequency region. The final
%  impedance characteristic is therefore constructed from three parts.
%  
%  If you need any additional information, please feel free to contact us at
%  andrej.debenjak@ijs.si
%  
%  Please cite as:
%  A. Debenjak, P. Boskoski, B. Musizza, J. Petrovcic, D. Juricic. Fast
%  measurement of proton exchange membrane fuel cell impedance based on 
%  pseudo-random binary sequence perturbation signals and continuous
%  wavelet transform. Journal of Power Sources (254), pp. 112-118, 2014. 
%% 
% Low frequency part from 0.05 to 5 Hz
disp('Low frequency part');
load meas_prbs__low.mat

nb = 20; %number of frequencies at which the CWT will be computed
cwt_data = zeros(6, nb); %prepare output data matrix

f_min = 0.05;  %lower frequency bound in Hz
f_max = 5;     %higher frequency bound in Hz
fs = 100000;   %sampling rate in Hz
decimate_factor = 10; %decimate factor for resempling (down-sampling)
f0 = 5;        %morlet central frequency

[Z f varA varPhi stat Wz c_u c_i ts] = ...
    fc_impedance_cwtft(u_ac, i_ac, f_min, f_max, fs, decimate_factor, f0, nb);

for kk=1:length(stat.real)
    cwt_data(1, kk) = stat.real(kk).median;
    cwt_data(2, kk) = stat.real(kk).prctile10;
    cwt_data(3, kk) = stat.real(kk).prctile90;

    cwt_data(4, kk) = stat.imag(kk).median;
    cwt_data(5, kk) = stat.imag(kk).prctile10;
    cwt_data(6, kk) = stat.imag(kk).prctile90;

    %cwt_data(7, kk) = stat.abs(kk).median;
    %cwt_data(8, kk) = stat.abs(kk).prctile10;
    %cwt_data(9, kk) = stat.abs(kk).prctile90;

    %cwt_data(10, kk) = stat.angle(kk).median;
    %cwt_data(11, kk) = stat.angle(kk).prctile10;
    %cwt_data(12, kk) = stat.angle(kk).prctile90;
end

Z_low = (cwt_data(1,:) + j*cwt_data(4,:));
f_low = f;

%% 
% Medium frequency part from 1 to 500 Hz
disp('Medium frequency part');
load meas_prbs__med.mat

nb = 40; %number of frequencies at which the CWT will be computed
cwt_data = zeros(6, nb); %prepare output data matrix

f_min = 1;     %lower frequency bound in Hz
f_max = 500;   %higher frequency bound in Hz
fs = 100000;   %sampling rate in Hz
decimate_factor = 1; %decimate factor for resempling
f0 = 15;       %morlet central frequency

[Z f varA varPhi stat Wz c_u c_i ts] = ...
    fc_impedance_cwtft(u_ac, i_ac, f_min, f_max, fs, decimate_factor, f0, nb);

for kk=1:length(stat.real)
    cwt_data(1, kk) = stat.real(kk).median;
    cwt_data(2, kk) = stat.real(kk).prctile10;
    cwt_data(3, kk) = stat.real(kk).prctile90;

    cwt_data(4, kk) = stat.imag(kk).median;
    cwt_data(5, kk) = stat.imag(kk).prctile10;
    cwt_data(6, kk) = stat.imag(kk).prctile90;

    %cwt_data(7, kk) = stat.abs(kk).median;
    %cwt_data(8, kk) = stat.abs(kk).prctile10;
    %cwt_data(9, kk) = stat.abs(kk).prctile90;

    %cwt_data(10, kk) = stat.angle(kk).median;
    %cwt_data(11, kk) = stat.angle(kk).prctile10;
    %cwt_data(12, kk) = stat.angle(kk).prctile90;
end

Z_med = (cwt_data(1,:) + j*cwt_data(4,:));
f_med = f;
     
%%
% High frequency part from 400 to 20000 Hz

disp('High frequency part');
load meas_prbs_high.mat

nb = 20; %number of frequencies at which the CWT will be computed
cwt_data = zeros(6, nb); %output data matrix

f_min = 400;      %lower frequency bound in Hz
f_max = 20000;    %higher frequency bound in Hz
fs = 200000;      %sampling rate in Hz
decimate_factor = 1; %decimate factor for resempling
f0 = 50;          %morlet central frequency

[Z f varA varPhi stat Wz c_u c_i ts] = ...
    fc_impedance_cwtft(u_ac, i_ac, f_min, f_max, fs, decimate_factor, f0, nb);

for kk=1:length(stat.real)
    cwt_data(1, kk) = stat.real(kk).median;
    cwt_data(2, kk) = stat.real(kk).prctile10;
    cwt_data(3, kk) = stat.real(kk).prctile90;

    cwt_data(4, kk) = stat.imag(kk).median;
    cwt_data(5, kk) = stat.imag(kk).prctile10;
    cwt_data(6, kk) = stat.imag(kk).prctile90;

    %cwt_data(7, kk) = stat.abs(kk).median;
    %cwt_data(8, kk) = stat.abs(kk).prctile10;
    %cwt_data(9, kk) = stat.abs(kk).prctile90;

    %cwt_data(10, kk) = stat.angle(kk).median;
    %cwt_data(11, kk) = stat.angle(kk).prctile10;
    %cwt_data(12, kk) = stat.angle(kk).prctile90;
end

Z_high = (cwt_data(1,:) + j*cwt_data(4,:));
f_high = f;
%%

figure();
plot(Z_low, 'b', 'linewidth', 2);
hold on
plot(Z_med, 'b', 'linewidth', 2);
plot(Z_high, 'b', 'linewidth', 2);
axis equal
grid on
title('Impedance characteristic [50 mHz -- 20 kHz]')
xlabel('Re [\Omega]')
ylabel('Im [\Omega]')


