clc
clear all
close all
Imax = 0.457
x = 0:(1/99):1;
x2 = 0:(3/299):3;
x2 = 0:(3/299/52e3*1e6):(3/52e3*1e6);
yintg_3 = 0:(3/299):3;
s = 0.2706;
s2 = 1-s;
Istr = 0:(1/99):1;
IC = 0:(1/99):1;
ICintg = 0:(1/99):1;
for i = 1:100
 Istr(i) = 0.1667;
end
krok = Imax/99/s2;
y = zeros(1,100);
y2 = Imax:-krok:0;
for i = 1:(s2*100)
 y(i) = y2(i);
end
%y(0:1:10) = 10
yintg = 0:(1/99):1;
%y(1:end) = 2;
y_3 = repmat(y,1,3);
for i = 1:99
 yintg(i+1) = yintg(i) + y(i)/99;
end
for i = 1:100
 IC(i) = y(i) - Istr(i);
end
for i = 1:99
 ICintg(i+1) = ICintg(i) + IC(i)/99/30e-6/52e3;
end
for i = 1:100
 ICintg(i) = ICintg(i) + 30;
end
y_3 = repmat(y,1,3);
yintg_3 = repmat(yintg,1,3);
IC_3 = repmat(IC,1,3);
ICintg_3 = repmat(ICintg,1,3);
Istr_3 = repmat(Istr,1,3);

plot (x2, y_3,"g");
hold on;
plot (x2, Istr_3,"k");
hold on;
 [hax, h1, h2] = plotyy (x2, IC_3, x2, ICintg_3);
legend ({"Proud sekundárním vinutím"
, "Střední proud odebíraný zátěží"
, "Proud tekoucí filtračním kondenzátorem"
, "Napětí na filtračním kondenzátoru"}
, "location", "southoutside");
 xlabel (hax(1), "t [\\mus]");
% title (hax(2), '');
 ylabel (hax(1), "I [A]");
 ylabel (hax(2), "U [V]");

%% THD

##%Following the proposed solution by Wayne, I will explain a method in order to %compute THD using the FFT.
##
##fs = 5200000 %6400; %sampling frequency
##Duree=1; %duration
##freq=50; %fundamental frequency of the signal
##t = linspace(0,Duree,Duree*fs);
##x = 100*sin(2*pi*50*t)+6*sin(2*pi*5*50*t)+5*sin(2*pi*7*50*t);
##
##%So we have three frequencies in the signal : 50Hz, 250 Hz and 350 Hz
##%The theoretical value of the THD is:
##
##THD_teorikoa = sqrt(6^2+5^2)/100
##% THD_teorikoa = 0.0781
##
##%We will use FFT method in order to compute this THD
##% The important element is to select a DFT_bin (DFT=Discrete Fourier Transform) %sufficiently accurate to be a multiple of the signal caracteristic frequencies (50, %250 and 350 Hz)
##
##DFT_bin = fs/length(ICintg); % = (sampling frequency)/(sample length)
##% DFT_bin = 1, then doing the FFT, the frquency 50 Hz, 350Hz and 350 Hz will be multiple of this DFT_bin frequency.
##% If you want more explanation you can put this words "Amplitude Estimation and %Zero Padding" in the research tool of matlab
##
##%Then now we use fft to compute the DFT and its power:
##m = length(ICintg); % Window length (number of samples)
##n = pow2(nextpow2(m)) ; % Transform length. This formula is used in order to accelerate the DFT computation, knowing that the computation is easiest if the %length is a multiple of two (prime number)
##Y = fft(ICintg,n)/n; %DFT computation
##f = fs/2*linspace(0,1,n/2+1); % Frequency range (the maximum of the frequency %range is nyquist frequency (Fs/2)+1
##
##% Plot single-sided amplitude spectrum.
##plot(f,abs(Y(1:n/2+1)))
##xlabel('Frequency (Hz)')
##ylabel('|Y|')
##title('{\bf Periodogram}')
##
##% Selection of the bin where the DFT is computed
##
##bin=fs/n; %The bin is computed taking into account that the frequency range is %divided in n bins. Then computation of the width of a bin will give us information, %where the interesting DFT at frequencies (50,250 and 350 Hz) are computed
##
##% Finally the computatio is:
##THD=sqrt(abs(Y(ceil((5*50)/bin+1)))^2+abs(Y(ceil((7*50)/bin+1)))^2)/sqrt(abs(Y(ceil(50/bin+1)))^2)
##
##% It's possible that an other formula will be better to compute the THD.
##hold on;
##yyaxis right
##plot (x2, ICintg_3);
## plot (x, y);
## hold on
## plot (x, yintg);
## hold on
## plot (x, IC);
## hold on
## plot (x, ICintg);
