%Fourier coefficients for BIPOLAR PWM 
% Copyright Daniel W. Hart, 2002
% See pp. 312-313 in the textbook Introduction to Power Electronics, D. Hart, 1997
% This m-file requires the function files alphacalc.m and alpdelta.m
clear
format compact
global m n b w ma

%%%%%%%%%%%% Enter values for these parameters %%%%%%%%%%%%%

mf=35           % freq. modulation ratio (odd integer)
ma=.9           % amplitude modulation ratio (1 or less)
vdc=100         % input dc source voltage
fo=60           % fundamental frequency (of sine reference)
nmax=mf*6 + 10  %number of terms to analyze in the Fourier series

r= 10    % resistance for series RL load
l=.02    % inductance in henries for series RL load

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


fc=fo*mf
w=2*pi*fo
p= (mf+1)/2;    % pulses per half cycle of sine

%calculate slope and intercept of each line segment of the triangular carrier signal
for n=1:2*p,
   m(n)=4*fc*(-1)^n ;
   b(n)=(2*n-2)*(-1)^(n+1);
end
%calculate alpha for each pulse

for n=1:p,
   alpha(n) = w*fzero('alphacalc',10e-3);
end

%calculate end of each pulse (alpha+delta)
for n=1:p,
   alpdelta(n)= w*fzero('alpdelta',1e-3);
end


%calc Fourier coefficients for Fbn

for n=1:nmax,
   fb(n)=0 ;
   for k=1:p,
      fb(n)=fb(n) + (2*vdc/(n*pi))*(cos(n*alpha(k)) - cos(n*alpdelta(k)));
   end
end

      for n=1:nmax,
         for k=1:(p-1) ;
            fb(n)=fb(n)+(2*vdc/(n*pi))*(cos(n*alpha(k+1)) - cos(n*alpdelta(k)));
         end
      end
      

for n=1:nmax,
   v(n)= abs(fb(n));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%LOAD CURRENT: SERIES R-L LOAD

%calculate current, impedance, power at each frequency
for n = 1:nmax
   irms(n) = (v(n)/sqrt(2)) / sqrt( r^2 + (n*w*l)^2) ;
   ipeak(n) = irms(n)*sqrt(2);
   power(n)=irms(n)^2*r;
   z(n) = sqrt(r^2+(n*w*l)^2);
end

I1rms = irms(1) %rms of fundamental

power_total=sum(power)

% calculate total harmonic distortion for the current
sumsquares=0;
for n = 2:nmax
   sumsquares = sumsquares + irms(n)^2;
end
sqrt(sumsquares);
thd = sqrt(sumsquares)/irms(1)
thdpercent = 100*thd

%% Graphs
subplot(2,1,1); bar(v), title(['voltage amplitude:   Vdc = ',num2str(vdc),', ma = ',num2str(ma),', mf = ',num2str(mf)])
subplot(2,1,2); bar(ipeak), title(['current amplitude:   THD = ',num2str(thdpercent,'%11.2g'),' %'] )