## preamp.m - calculations on a two stage preamp for a multi-turn,
## air-core solenoid loop antenna for the reception of signals below
## 30kHz.
##
## Usage: preamp(verbose)
##   if verbose <> 0 be noisey
##   default: don't be noisey
##
function [TFpreamp]=preamp(verbose=0)
##
disp(verbose)
disp("");
disp("VLF Pre-Amplifier Design");
disp("This example covers the design of a pre-amplifier for use in receiving");
disp("radio frequencies below 30kHz.");
disp("See http://www.vlf.it for details of Narural Radio Sources");
disp("");
disp("The preamp consists of two AD797 op amps and a low pass filter. With");
disp("biasing and blocking capacitors omitted, three blocks remain.");
disp("");
disp("");
disp("          Gain = 10");
disp("        +-------------+");
disp("        |             |      -- Low Pass Filter --");
disp("     ---+ p           |");
disp(" Loop   |   Stage 1   +--+----R3--+--R4--+--R5--+---> To Stage 2");
disp("        |   Amplifier |  |        |      |      |");
disp("     -+-+ n           |  |        C1     C2     C3");
disp("      | |             |  |        |      |      |");
disp("      | +-------------+  |       Gnd    Gnd    Gnd");
disp("      |                  |");
disp("      +----+---R2--------+");
disp("           |");
disp("           R1");
disp("           |");
disp("          Gnd");
disp("");
disp("");
disp("");
disp("                   Gain = 10");
disp("                 +-------------+");
disp("                 |             |");
disp("            Gnd--+ p           |");
disp("                 |   Stage 2   +--+----R8--+----> Output");
disp("                 |   Amplifier |  |        |");
disp("      From >---+-+ n           |  |        R9");
disp("     Filter    | |             |  |        |");
disp("               | +-------------+  |       Gnd");
disp("               |                  |");
disp("               +----+---R6--------+");
disp("                    |");
disp("                    R7");
disp("                    |");
disp("                   Gnd");
disp("");
disp("");
disp("R1 and R2 profide feedback to control the gain of Stage 1.");
disp("R3 through R5 with C1 through C2 form a low pass filter to limit the bandwidth.");
disp("R6 and R7 profide feedback to control the gain of Stage 1.");
disp("R8 and R9 provide impedance matching to the cable and/or receiver,");
disp("possibly a PC sound card.");
disp("");
##
ans=input("Press Return to continue:");
##
disp("");
disp("");
disp("The graphs in the ad797 datasheet reveal the following parameters:");
disp("  Gain, a0 = 1e7 ");
disp("  Pole 55 Hz");
disp("  Pole 1 MHz");
disp("  Zero 4.3 MHz");
disp("");

a0 = show("a0 = 1e7");
p1 = show("p1 = 55");
p2 = show("p2 = 1e6");
z1 = show("z1 = 4.3e6");

disp("");
disp("The open loop transfer function of an op amp with m zeros and n");
disp("poles is expressed in the form:");
disp("  tf = dcgain * zero_expressions / pole_expressions");
disp("where ");
disp("  zero_expressions = (1+s/z1) * (1+s/z2) * ... * (1+s/zm) ");
disp("  pole_expressions = (1+s/p1) * (1+s/p2) * ... * (1+s/pn)");
disp("where ");
disp("  z1 ... zm are the m zeros");
disp("  p1 ... pn are the n poles");
disp("");
disp("or, in the case of 3 zeros and 4 poles:");

s = show("s = tf('s')")
disp("");
TFopen = show("TFopen = a0 * (1+s/z1) / (1+s/p1) / (1+s/p2)")
disp("");
TFopen_norm = show("TFopen_norm = minreal(TFopen)")

##
ans=input("Press Return to continue:");
##

disp("");
disp("Note: The difference between the op amp expression and the usual");
disp("Zero-Pole-Gain expression is in the modification of the dcgain");
disp("parameter.  The dc gain argument to zpk() is modified by the zeros and poles,");
disp("");

Azpk = show("Azpk = zpk([-z1],[-p1,-p2],1e7*p1*p2/z1)")

##
ans=input("Press Return to continue:");
##
disp("");
disp("The bode plot of these two developments of the open loop transfer function");
disp("produce identical results.  Ane, the plots show the same shape as the");
disp("graphs in the datasheet.");
disp("");

show("figure 1");
show("bode(TFopen)");
show("subplot(2,1,1)");
show("title('Equation Bode Diagram')");
show("figure 2");
show("bode(Azpk)");
show("subplot(2,1,1)");
show("title('ZPK Bode Diagram')");

disp("");
disp("Two Bode Diagrams should be visible, possibly ovelaid.");
disp("");
##
ans=input("Press Return to close the plots and continue:");
##

close all;
disp("");
disp("The normalized step response of the ad797 is:");

disp("");
TFnorm = show("TFnorm = TFopen/dcgain(TFopen)")
disp("");
show("step(TFnorm,'b')");
show("title('AD797 Normalized Open-Loop Step Response')");
show("ylabel('Normalized Amplitude')");

##
ans=input("Press Return to close the plot and continue:");
##
close all;
disp("Design Stage 1 of the preamp.");
disp("Resistors R1 and R2 form a feedback system to control the gain of Stage 1");
disp("This feedback system returns a portion of the output to the negative input.");
disp("This is normally expressed as:");
disp("  Vfb = Vout * R1 / (R1 + R2)");
disp("So, the transfer function of the feedback network is:");
disp("  tf = Vfb / Vout = R1 / (R1 + R2)");
disp("The effects of the AD797 gain on the input and the feedback may be represented");
disp("as TFstage1 = Vout/Vp = dcgain / (1 + dcgain * TFfeedback).");
disp("If dcgain is sufficiently large, this reduces to TFstage1 = 1 / TFfeedback.");
disp("The dcgain of the AD797 is >> 1, so, the feedback completely controls the");
disp("output and variations in the dcgain will not effect the Stage gain.");
disp("");
disp("The feedback is added to the AD797 using the feedback function");

Gfb = show("Gfb = 10");
b   = show("b = 1 / Gfb");
R1  = show("R1 = 10e3");
R2  = show("R2 = R1 * (1/b - 1)")
disp("");
TFstage1 = show("TFstage1 = feedback(TFopen,b)")
disp("");

show("bodemag(TFopen,'r',TFstage1,'b')");
show("legend('Open Loop Gain (TFopen)','Closed Loop Gain (TFstage1)')");

disp("The use of negative feedback to reduce the low-frequency (LF) gain");
disp("has led to a corresponding increase in the system bandwidth (defined");
disp("as the frequency where the gain drops 3dB below its maximum value).");
disp("");
disp("With this feedback, we have a gain of 10, or 20db up to 10MHz,");
disp("far more than the frequency range of interest.");
disp("");
##
ans=input("Press Return to close the plot and continue:");
##
close all;
disp("Since the gain is now dominated by the feedback network, a useful");
disp("relationship to consider is the sensitivity of this gain to variation");
disp("in the op amp's natural (open-loop) gain.");
disp("");
disp("Before deriving the system sensitivity, however, it is useful to");
disp("define the loop gain, L(s)=a(s)b(s), which is the total gain a signal");
disp("experiences traveling around the loop:");
disp("");
disp("Sensitivity = partial(TFstage1/TFopen)*TFopen/TFstage1");
disp("or S(s) = 1 / (1 + TFopen(s) * TFstage1(s))");
disp("or S(s) = 1 / (1 + L(s)), which has the same form as feedback");
disp("So, use the feedback function to develop the sensitivity.");
disp("");

L = show("L = TFopen * b")
disp("");
S = show("S = feedback(1, L)")
show("figure 1");
show("bodemag(TFstage1,'b',S,'g')");
disp("");

disp("The very small low-frequency sensitivity (more than -100 dB) indicates a");
disp("design whose closed-loop gain suffers minimally from open-loop gain");
disp("variation. Such variation in a(s) is common due to manufacturing");
disp("variability, temperature change, etc.");
disp("");
disp("You can check the step response of A(s) using the STEP command:");

show("figure 2");
show("step(TFstage1)");

disp("");
disp("The stability margin can be analyzed by plotting the loop gain, L(s) with");
disp("the margin function.");
disp(" ");
disp("This plot may display some warning messages, you can safely ignore them.");
disp(" "); fflush(stdout);
show("margin(L)");
disp(" "); fflush(stdout); fflush(stderr);

##
ans=input("Press Return to close the plot and continue:");
##

disp("");
disp("The plot indicates a phase margin of less than 3 degrees.  Stage 1 needs");
disp("to be compensated to increase this to an acceptible level,");
disp("more than 45 degrees, if possible.");
disp("");
disp("Feedback Lead Compensation");
disp("A commonly used method of compensation in this type of circuit is");
disp("feedback lead compensation. This technique modifies b(s) by adding");
disp("a capacitor, C, in parallel with the feedback resistor, R2.");
disp("The capacitor value is chosen so as to introduce a phase lead to b(s)");
disp("near the crossover frequency, thus increasing the amplifier's phase");
disp("margin.");
disp("The new feedback transfer function is shown below.");
disp("You can approximate a value for C by placing the zero of b(s) at the");
disp("0dB crossover frequency of L(s):");
disp("");

[Gm,Pm,Wcg,Wcp] = show("[Gm,Pm,Wcg,Wcp] = margin(L)");
C = show("C = 1/(R2*Wcp)")
show("if (C < 1e-12) disp('The calculated value of C is very small.  Look at a range of values.'); endif");

disp(" ");
disp("The next plots take some time...");
##
ans=input("Press Return to close the plot and continue:");
##
close all;
disp(" ");

K = show("K = R1/(R1+R2)");
C = show("C = [10:10:200]*1e-12");
b_array = show("b_array = arrayfun (@(C) tf([K*R2*C, K], [K*R2*C, 1]), C,'uniformoutput',false)");
A_array = show("A_array = cellfun(@feedback,{TFopen},b_array, 'uniformoutput', false)");
L_array = show("L_array = cellfun(@sys_prod,{TFopen},b_array, 'uniformoutput', false)");
S_array = show("S_array = cellfun(@feedback,{tf(1,1)},L_array, 'uniformoutput', false)");
disp(" "); fflush(stdout);
[Gm,Pm,Wcg,Wcp] = show("[Gm,Pm,Wcg,Wcp] = cellfun(@margin,L_array)");
disp(" ");

close all
show("figure 1");
show("step(TFstage1,'r',A_array{:})");
show("figure 2");
show("bode(TFstage1,A_array{:})");
show("figure 3");
show("plot(C,Pm)");
show("grid");
show("xlabel('Compensation Capacitor, C (pF)')");
show("ylabel('Phase Margin (deg)')");
show("figure 4");
show("step(A_array{C==50e-12},'r',A_array{C==100e-12},'b',A_array{C==200e-12},'g')");
show("legend('Compensated (4 pF)','Compensated (10 pF)','Compensated (20 pF)')");

disp(" ");
##
ans=input("Press Return to close the plot and continue:");
##
close all;
disp("");
disp("");
disp("          Gain = 10");
disp("        +-------------+");
disp("        |             |");
disp("     ---+ p           |");
disp(" Loop   |   Stage 1   +--+---->");
disp("        |   Amplifier |  |");
disp("     -+-+ n           |  |");
disp("      | |             |  |");
disp("      | +-------------+  |");
disp("      |                  |");
disp("      +----+------R2-----+");
disp("           |             |");
disp("           +-----Ccomp---+");
disp("           |");
disp("           |");
disp("           R1");
disp("           |");
disp("          Gnd");
disp("");
disp("The selected compensation is 10pf.");

TFcomp = show("TFcomp = A_array{C==100e-12}");
show("bode(TFopen,'b',TFstage1,'g',TFcomp,'r')");
show("legend('TFopen','TFstage1','TFcomp')");

disp("");
##
ans=input("Press Return to close the plot and continue:");
##
close all;
disp("Low Pass Filter Design");
disp("");
disp("The low pass filter is composed of three equal sections.");
disp("Develop one section and put three in series.");

C = show("C = 20e-9");
R = show("R = 1000");
TFsection = show("TFsection=tf( [1] , [C*R,1] )");
TFfilter = show("TFfilter=series( series( TFsection , TFsection ) , TFsection)");
if (verbose)
  TFfilter
endif;

disp("");

TFpreamp = show("TFpreamp = series(series(TFcomp,TFfilter),TFcomp)");
if (verbose)
  TFpreamp
endif;
show("figure 1");
show("bode(TFpreamp,{1,1e5})");
show("figure 2");
show("margin(TFpreamp)");

disp("");
disp("Two plots are displayed, possibly overlaid.");
##
ans=input("Press Return to close the plot and continue:");
##
disp("As can be seen from the plots, the gain margin is almost 30db.");
disp("The phase margin is 230 degrees.");
disp("");
disp("Use 'close all' to close the plots.");
##
endfunction;

function [x]=sys_prod(S1,S2)
  x=S1*S2;
endfunction;

function [r,s,t,u] = show(str)
  disp (str);
  switch nargout
  case 1
    r = evalin ("caller", str);
  case 2
    [r,s] = evalin ("caller", str);
  case 3
    [r,s,t] = evalin ("caller", str);
  case 4
    [r,s,t,u] = evalin ("caller", str);
  otherwise
    evalin ("caller", str);
  endswitch;
endfunction;