%This program computes data streams in a PWM power amp with feedback 
clear  % this is necessary to clear variables from prior runs
%***************************************************************
%***************************************************************
%*** Below is the section for the user to enter and edit data. 
%***************************************************************
%***************************************************************
N=16384 ;% number of time points in simulation, hence # of frequency points
% N needs to be a power of 2 for the simulation to run fast. 
% The student version of matlab limits you to 16238 (2^14) time points. 
% much better accuracy will be obtained if more than 16238 time points are used
Vp= 10 ; % *********this is the peak (not peak-peak) output voltage of the switching stage
beta=0.1 ; % *********this is the feedback factor R2/(R1+R2).
% From feedback theory, the closed-loop gain will be very close to 1/beta
% so the maximum input voltage you should apply should be a little less than Vp*beta.
f_tri=500E3 %triangle wave frequency in Hz
v_tri_mag=1 ;  % magnitude of triangle wavewaveform in volts
f_clock=f_tri*8; %clock frequency in Hz. Must be K*f_tri, where K is an integer.
% fclock sets how often the simulation makes a time step
% f_clock is NOT the triangle wave frequency. fclock is necessary becuase the program must
% make calculations at a finite number of steps in time.  So note that the in this simulation 
% that the pulse width modulation is now fully continuous in this simulation. If you make N larger
% than 16384, you can make f_clock larger. If you make fclock larger *without* making N larger,
% there will not be enough time points in one cycle of the signal, and the FFT will be noisy. 
loop_bandwidth=100E3% this is the feedback loop bandwidth
% the feedback loop bandwidth should be as high as possible for low distortion,
% but at MOST 1/5 of the triangle wave frequency. Problems arise otherwise
f_unity=loop_bandwidth*v_tri_mag/Vp/beta% unity gain frequency of the 2 integrators
% this relationship arises because the voltage gain of the PWM stage is (Vp/v_tri_mag)
%
f_zero=1E3 ; % frequency of the zero in the second integrator
% *********this should be less than the feedback loop bandwidth by at least
% a factor of 2, and probably by a factor of 10. 
f_cutoff=50E3 ; % cutoff frequency of the output filter
% usually set to somewhat higher than the maximum input signal frequency
n_signal_cycles=16; 
%  the signal frequency is given by f_signal=n_signal_cycles/tfinal
%  specifying the signal frequency in this fashion forces an interger # of
%  signal periods in the fft, which is necessary to obtain a clean fft. 
%  So, if you change N or f_tri, you will need to change n_signal_cycles in order
%  to get the signal frequency you want. 
v_in = 0.8 ;% the input signal voltage amplitude, peak, not RMS
gm1=2E-3 ; %transconductance of the first Gm element
gm2=2E-3 ; %transconductance of the second Gm element
c1=gm1*0.159/f_unity  % first integration capacitor
Rzero=1/gm2 ; % the resistor which introduces a zero in the second integrator
% this forces the 2nd integrator gain to unity at frequencies *above* the zero
c2=0.159/f_zero/Rzero  % second integration capacitor
Idac1_val=beta*gm1*Vp ; % 
% If the feedback is in the form of a voltage divider of
% attenuation beta, feeding back to the *input* of the first gm stage,
% then Idac1_val should be given by Idac1_val=beta*gm1*Vp
% This is the configuration most probably used by 137b students
%***************************************************************
%***************************************************************
%***************************************************************
%***************************************************************
%The program now sets up time and frequency increments for the simulation. 
f_nyquist=2*f_clock ;
fhigh=2*f_nyquist ;
tfinal=1/fhigh*(N-1)^2/N ;% 
f_signal=n_signal_cycles/tfinal % 
timestep=tfinal/(N-1); % timestep
f=linspace(0,fhigh,N);
delta_f=1/tfinal  % frequency step
t=linspace(0,tfinal,N);
N_clock=f_clock*tfinal ;  % number of data bits in simulation
t_clock=1/f_clock ; % the clock period. 
%*************************************************************************
%*************************************************************************
% the values of these variables are first set to very small values
% before running the finite-time simulations. 
%Seems to be necessary to avoid a progam bug.
% The real values are subsequently calculated in the next block down
v1=1E-9*sin(2*pi*f_clock*t); % v1 is the voltage on the first integration capacitor C1 
v2=1E-9*sin(2*pi*f_clock*t); % v2 is the voltage on the second integration capacitor C2
v3=1E-9*sin(2*pi*f_clock*t); % v3 is the voltage at the output of the second integrator 
v_tri=1E-9*sin(2*pi*f_clock*t); % v_clock is the triange wave generator
v_data=1E-9*sin(2*pi*f_clock*t); 
v_data2=1E-9*sin(2*pi*f_clock*t); 
%v_data(n_count)=1 ; % v_data is the output of the latched  comparator
Idac1=1E-9*sin(2*pi*f_clock*t); % the output of the first DAC
Idac2=1E-9*sin(2*pi*f_clock*t); % the output of the second DAC
%*************************************************************************
%*************************************************************************
% now we simulate the feedback loop in finite time steps
% this is the core of the simulation
% First apply an input signal: a sinewave  
input_signal=v_in*sin(2*pi*f_signal*t); % 
v_tri=v_tri_mag*( sin(2*pi*f_tri*t)+sin(4*pi*f_tri*t) ); % 
time=0 ;
n_count=1; % intializing the count variable
time_int=0; % this is accumulated time within a clock cycle...see below
% we now have a if/else/for loop which increments time in small steps
for n_count=2:N ;
   squarewave=1 ;
   n_countless=n_count-1;
   time=n_count*timestep; % increment time by one timestep
   input_signal(n_count)=v_in*sin(2*pi*f_signal*time);% the input signal
   v_tri(n_count)=cos(2*pi*f_tri*time); % the triangle wave: fundamental
   v_tri(n_count)=v_tri(n_count)+(1/9)*cos(3*2*pi*f_tri*time); % the triangle wave: adding the 3rd harmonic
   v_tri(n_count)=v_tri(n_count)+(1/25)*cos(5*2*pi*f_tri*time); % the triangle wave: adding the 5th harmonic
   v_tri(n_count)=v_tri(n_count)+(1/49)*cos(7*2*pi*f_tri*time); % the triangle wave: adding the 7th harmonic
   v_tri(n_count)=v_tri(n_count)+(1/81)*cos(9*2*pi*f_tri*time); % the triangle wave: adding the 9th harmonic
   v_tri(n_count)=v_tri(n_count)+(1/121)*cos(11*2*pi*f_tri*time); % the triangle wave: adding the 11th harmonic
   v_tri(n_count)=v_tri(n_count)+(1/169)*cos(13*2*pi*f_tri*time); % the triangle wave: adding the 13th harmonic
   v_tri(n_count)=v_tri(n_count)+(1/225)*cos(15*2*pi*f_tri*time); % the triangle wave: adding the 15th harmonic
   v_tri(n_count)=v_tri(n_count)+(1/289)*cos(17*2*pi*f_tri*time); % the triangle wave: adding the 17th harmonic
   v_tri(n_count)=v_tri_mag*v_tri(n_count)/1.2; % the triangle wave: scaling for the correct amplitude
   if v_data(n_countless)> 0 ; % setting feedback signal according to the comparator output
      Idac1(n_count)=Idac1_val;
   else
      Idac1(n_count)=-1*Idac1_val;
   end
   % summing currents at C1, and then doing I = c dV/dt to find the voltage on c1
   I_nodeone(n_count)= gm1*input_signal(n_count)-Idac1(n_count);
   n_count_less=n_count-1;
   v1(n_count)=v1(n_count_less)+I_nodeone(n_count)*timestep/c1;
   % finding the current on C2 and then doing I = c dV/dt to find the voltage on c2
   I_nodetwo(n_count)= gm2*v1(n_count);
   v2(n_count)=v2(n_count_less)+I_nodetwo(n_count)*timestep/c2;
   % adding the voltage drop on the resistor Rzero 
   v3(n_count)=v2(n_count)+Rzero*gm2*v1(n_count) ; 
   time_int=time_int+timestep; % time accumulated within the clock period
   if time_int < t_clock ;% if the accumulated time within a bit is less than a clock period
      v_data(n_count)=v_data(n_count_less) ;% then the data keeps its same value
      v_data2(n_count)=squarewave*v_data(n_count);
   else 
      % otherwise the data changes according to the new comparator input
      v_data(n_count)=Vp*sign(v3(n_count)-v_tri(n_countless));%this line for 2 integrators
      %v_data(n_count)=Vp*sign(v1(n_count)-v_tri(n_countless));% this line for 1 integrator
      time_int = time_int -t_clock ; % and the accumlated time is re-incremented (new clock!)
   end
end
%*************************************************************************
%*************************************************************************
% now we FFT the data stream, and then filter it so that we can see the 
% resulting low-pass-filtered output signal 
%
% first we must shuffle the frequencies in the modulation filter
% to get them in the strange order the FFT requires
alpha=0.51 ;
n_count=1;
for n_count = 1:N
   if n_count > N/2
      f_c=1/tfinal*(N-1)/N*(n_count-N);
   else
      f_c=1/tfinal*(N-1)/N*(n_count-1);
   end
   H_filter(n_count)=raised_cos_baseband_filter2(f_c,f_cutoff,alpha); % root-raised-cos filter
   % the filter I have used for convenience is a mathematically idealized raised-cosine filter.
   % your real circuit has the low pass filter determined by the values of L and C you use
   % on the output. 
end
%
%*************************************************************************
%now we can filter the modulated signal

v_data_f=fft(v_data) ;% first FFT the data stream
input_signal_f=fft(input_signal); % also handy to have the fft of the input signal 
%
%
n_count=1;
for n_count=1:N % then multiply it with the filter transfer function. 
   filt_data_f(n_count)=v_data_f(n_count)*H_filter(n_count) ;  
end
filt_data_t=ifft(filt_data_f); % and then IFFT to get the time waveform
%
%
in_spec = fft(input_signal.*hann(N));
out_spec = fft(v_data.*hann(N));
%
%
%*************************************************************************
%*************************************************************************
%*************************************************************************
%*************************************************************************
%*************************************************************************
%
%Plotting calls are below
%
close all
%***************************************************************
% **************** first plot window ***************************
scrsz = get(0,'ScreenSize');
%figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
% left bottom width height
fg1=figure('Position',[0 0 scrsz(3)/1.1 scrsz(4)/1.1]);
%
% this plots the first integrator output voltage
ndatasplot=30; % enter here the number of signal cyles you want plotting
t_sigplot_stop=ndatasplot*t_clock;
plot(t,v1)
xlabel('time, sec')
%axis([0 t_sigplot_stop -2.5 2.5])
%axis([0 tfinal -2 2])
title('first integrator output voltage: MAKE SURE IT IS LESS THAN THE CIRCUIT Voltage CLIPPING LIMITS')
%
%***************************************************************
% **************** 2nd plot window ***************************
scrsz = get(0,'ScreenSize');
%figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
% left bottom width height
fg2=figure('Position',[0 0 scrsz(3)/1.1 scrsz(4)/1.1]);
% this plots the output voltage of the 2nd integrator
ndatasplot=30; % enter here the number of signal cyles you want plotting
t_sigplot_stop=ndatasplot*t_clock;
plot(t,v3)
xlabel('time, sec')
%axis([0 t_sigplot_stop -2.5 2.5])
%axis([0 tfinal -2 2])
title('second integrator output voltage: MAKE SURE IT IS LESS THAN THE CIRCUIT Voltage CLIPPING LIMITS')
%
%***************************************************************
% **************** 3rdplot window ***************************
scrsz = get(0,'ScreenSize');
%figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
% left bottom width height
fg3=figure('Position',[0 0 scrsz(3)/1.1 scrsz(4)/1.1]);
%
% this plots the ouput voltage unfiltered and filtered
ndatasplot=300; % number of data cyles you want plotting
t_sigplot_stop=ndatasplot*t_clock;
plot(t,v_data,t,filt_data_t)
xlabel('time, sec')
axis([0 t_sigplot_stop -1.1*Vp 1.1*Vp])
%axis([0 tfinal -2 2])
title('ouput voltage, unfiltered digtal swiching pattern and filtered analog output voltage')
%
%***************************************************************
% **************** 4th plot window ***************************
scrsz = get(0,'ScreenSize');
%figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
% left bottom width height
fg4=figure('Position',[0 0 scrsz(3)/1.1 scrsz(4)/1.1]);
%
% this plots input signal voltage
ndatasplot=30; % enter here the number of signal cyles you want plotting
t_sigplot_stop=ndatasplot*t_clock;
plot(t,input_signal)
xlabel('time, sec')
%axis([0 t_sigplot_stop -2.5 2.5])
%axis([0 tfinal -2 2])
title('input signal voltage ')
%
%
%***************************************************************
%***************************************************************
% **************** 5th plot window ***************************
fg5=figure('Position',[0 0 scrsz(3)/1.1 scrsz(4)/1.1]);
%********************************************
% this plots the spectrum of the data
%semilogx( f ,20*log10(2*abs(v_data_f)/N+1e-90),f, ...
%   20*log10(2*abs(input_signal_f)/N+1e-90));
semilogx( f ,20*log10(2*abs(in_spec)/N+1e-90),f, ...
   20*log10(2*abs(out_spec)/N+1e-90));
%axis([ 0 5*f_clock -50 -10 ])
%axis([ 0 log10(fhigh/2) -100 -10 ]) 
axis([ 0 fhigh/2 -100 20*log10(abs(Vp)) ]) 
%set('XScale', 'log)
xlabel('Frequency, Hz')
ylabel('Power, dB_V , with  delta-f Hz FFT frequency resolution ')

% you need to work out the FFT frequency resolution !!!!!
title('input and output signals in the frequency domain') 
%
%***************************************************************
%***************************************************************
% **************** 6th plot window ***************************
scrsz = get(0,'ScreenSize');
%figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
% left bottom width height
fg6=figure('Position',[0 0 scrsz(3)/1.1 scrsz(4)/1.1]);
%
% 
ndatasplot=10; % enter here the number of signal cyles you want plotting
t_sigplot_stop=ndatasplot*t_clock;
plot(t,I_nodeone)
xlabel('time, sec')
%axis([0 t_sigplot_stop -2.5 2.5])
%axis([0 tfinal -2 2])
title('first gm stage output current: MAKE SURE IT IS LESS THAN THE CIRCUIT CURRENT CLIPPING LIMITS')
%***************************************************************
% **************** 7th plot window ***************************
scrsz = get(0,'ScreenSize');
%figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
% left bottom width height
fg7=figure('Position',[0 0 scrsz(3)/1.1 scrsz(4)/1.1]);
%
% this plots second gm stage output current
ndatasplot=10; % enter here the number of signal cyles you want plotting
t_sigplot_stop=ndatasplot*t_clock;
plot(t,I_nodetwo)
xlabel('time, sec')
%axis([0 t_sigplot_stop -2.5 2.5])
%axis([0 tfinal -2 2])
title('second gm stage output current: MAKE SURE IT IS LESS THAN THE CIRCUIT CURRENT CLIPPING LIMITS')
%
%***************************************************************
% **************** 8th plot window ***************************
scrsz = get(0,'ScreenSize');
%figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
% left bottom width height
fg8=figure('Position',[0 0 scrsz(3)/1.1 scrsz(4)/1.1]);
%
% this plots second gm stage output current
ndatasplot=150; % enter here the number of signal cyles you want plotting
t_sigplot_stop=ndatasplot*t_clock;
plot(t,v_tri,t,v1,t,v_data/Vp )
xlabel('time, sec')
axis([0 t_sigplot_stop -2.5 2.5])
%axis([0 tfinal -2 2])
title('triangle wave (blue), comparator input (green), normalized comparator output (red)')
%