%  This script runs an example of a load line analysis for a MOSFET amp
%  to:
%   1. find the operating point for Vgs through iteration
%   2. then compute a table of input (Vgs) and output (Vds)amplitude values
%      to get the transfer characteristic (TC) curve for the amp.
%   3. then apply the resulting TC curve to an input sinusoid with
%      increasing amplitude to illustrate distorion.  The sine wave
%      will be plotted and played in the demonstration.
%  The functions nmos.m and ampdist.m are needed to run this script
%
%   Written by Kevin D. Donohue (donohue@engr.uky.edu) 10/2003

%  Set Parameters:  Operating point will be set to half VDD
K=.5; vto = 1.8;  %  Nmos parameters
W=1; L=1; KP=2*K; %  Nmos parameters
VDD=15;  RS=400; RD=1e3;  %  Load line parameters
idsmax = VDD/(RD+RS);    % Maximum Load line value on Drain current axis
                         %  Operation point will be the max drain current divided by 2.
vds = [0:.05:VDD];  %  Create X-Axis

idsll = -vds/(RD+RS) + VDD/(RD+RS);  %  Generate Load Line

err = .2e-3;   %Set initial error for quiesent point to get while loop started
toler = .05e-3;  %  Set tolerance value for stopping rule on while loop
incgs = .01;  %  Set increment for vgs to find intersection with load line.
vgs = vto;         %  Initialize vgs to threshold voltage so initial guess will below desired operating point
%  Set flag to denote when guess goes above the the operating point in
%  order to reduce interation interval.
passflag = 0;  %  If zero implies the last step is below the desired value, 1 implies last step was above

%  Set while loop to run until error is below tolerance
while err >= toler
	ids = nmos(vds,vgs,KP,W,L,vto);  % Compute characteristic curve
	[edum, inderr] = min(abs(idsll - ids));  %  Find intersection with loadline
    err =abs(idsmax/2 - idsll(inderr(1)));  %  Find error between desired current and actual
    % Check to see if we are above or below the target and make adjustments to move closer to desired operating point
	if idsmax/2>ids(inderr)  %  If below target value ....
        vgs = vgs + incgs   %  Still below threshold so increase
        if passflag == 1
            incgs = incgs/2;   %  If we just came from above the threshod cut increment in half for more resolution
            passflag = 0;      %  reset flag
        end
    else                   %  If above target value
        vgs = vgs - incgs   % Still above threshold so decrease 
        if passflag == 0
            incgs = incgs/2;  % If we just came from below the threshold cut increment in half for more resolution
            passflag = 1;  % reset flag
        end
	end
%   Comment out the next 4 lines to stop the while loop from being interupted by plots and pauses
    figure(1)
    plot(vds, idsmax*ones(size(vds))/2, 'c--', vds,idsll,'k:',vds,ids,'r')  %  Check plot along the way
    disp(['Iteration in progress'])
    pause(.1)
   end
vgsq = vgs;  %  Set quiesent Vgs to last result of interation
disp([' The operating point for Vgs is ' num2str(vgsq)])

%  Compute corresponding output quiesent voltage amplitude at output
ids = nmos(vds,vgsq,KP,W,L,vto);  % Compute transfer characteristic (TC) curve at quiesent
[err, inderr] = min(abs(ids - idsll));   %  Find intersection with load line and TC
vdsq = vds(inderr(1));  %  Output voltage quiesent
disp([' The operating point for Vds is ' num2str(vdsq)])
idsq = ids(inderr(1));
disp([' The operating point for Ids is ' num2str(idsq)])


disp([' Now compute transfer charactersitic between amp input and output '])

%   Now compute array for maping the input to the output of the amplifier

inarray = [0:.001:2*vgsq];         %   AC input array amplitude sweep
                                   %   range should include quiesent Vgs
%  Loop to compute each point on intersection of load line for input
%  vortage amplitude sweep for output VDS
outarray = zeros(1,length(inarray)); %  Initalize output array
for karry = 1:length(inarray)
         vgsdcac= inarray(karry);   %  Vgs level with AC and DC energy
         ids = nmos(vds,vgsdcac,KP,W,L,vto);  % Compute characteristic curve for that vgs
         [err, inderr] = min(abs(ids - idsll));   %  Find closest point between load line and TC
         outarray(karry) = vds(inderr(1));      %  Assign VDS as output array correspond to that Vgs value
end
%  Subtract quiesent Vgs offset from input array to result in AC input only component
insigac = inarray - vgsq;
%  Subtract quiesent Vds offset from input array to result in AC output only component
outsigac = outarray - vdsq; 
%  Plot transfer characteristics for AC voltage gain
figure(2)
plot(insigac,outsigac)
xlabel('Input AC voltage Amplitude')
ylabel('Output AC voltage Amplitude')
title('AC transfer characteristic of amplifier')
disp(['Hit any key to continue to hear examples of sounds played through amplifier'])
pause
%  Now test the input and output relationship of the amplifier:

%  Creat a unit 300 Hz sine wave sampled at 8000 Hz and pass it through the amp

fs = 8000;  %  Sampling frequency
t= [0:3*fs-1]/fs;    %  Create a time axis for signal for 3 seconds
sigin = sin(2*pi*t*300);   %  Create unit sine wave
a = [.01, .05, .15];  %Set up amplitude scales for input voltage
%  Loop to distort and play sound
for k=1:3
    sigout = ampdist(sigin*a(k),insigac,outsigac);  %  Distort sound by mapping amplitudes through TC
    %  Plot original and amplified signal on same scale to observe distortion
    figure(2+k); plot(t(1:100),sigin(1:100)/max(abs(sigin)),'r',t(1:100),-sigout(1:100)/max(abs(sigout)),'b')
    title([' Compare Scaled Input (red) at amplitude ' num2str(a(k)), ' output (blue) for Distortion'])
    %  Play both sounds consecutively
    %  Concatenate scale input and output signals in one vector 
    wavsig = [sigin/(max(abs(sigin))+eps), sigout/(max(abs(sigout))+eps)];
    soundsc(wavsig,fs);
    %  Pause for user key press to go on to next sound if not at the end
    if k~=3
        disp(['Hit any key to continue to next sound'])   
        pause
    end
end