## preamp.m - calculations on a two stage preamp for a multi-turn,
## air-core. solenoid loop antenna for the reception of signals below
## 30kHz.
function [TFpreamp]=preamp()
##1; ## not a function file
##
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 = 1e7;
p1 = 55;
p2 = 1e6;
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 = tf('s');
TFopen = a0 * (1+s/z1) / (1+s/p1) / (1+s/p2);
disp("");
TFopen_norm = minreal(TFopen)
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("");
disp("zpk([-z1],[-p1,-p2],1e7*p1*p2/z1)");
disp("");
Azpk=zpk([-z1],[-p1,-p2],1e7*p1*p2/z1)
##
ans=input("Press Return to continue:");
##
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("");
figure 1;
bode(TFopen);
subplot(2,1,1);
title("Equation Bode Diagram");
figure 2;
bode(Azpk);
clear eqn;
subplot(2,1,1);
title("ZPK Bode Diagram");
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:");
TFnorm=TFopen/dcgain(TFopen);
step(TFnorm,'b');
title('AD797 Normalized Open-Loop Step Response');
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 = 10;
b = 1/Gfb;
R1 = 10e3;
R2 = R1 * (1/b - 1)
TFstage1 = feedback(TFopen,b);
disp("");
bodemag(TFopen,'r',TFstage1,'b')
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 = TFopen * b;
S = feedback(1, L)
figure 1;
bodemag(TFstage1,'b',S,'g');
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:");
figure 2;
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.");
##
margin(L);
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] = margin(L);
C = 1/(R2*Wcp)
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 = R1/(R1+R2);
C = [10:10:200]*1e-12;
b_array=arrayfun (@(C) tf([K*R2*C, K], [K*R2*C, 1]), C,'uniformoutput',false);
A_array=cellfun(@feedback,{TFopen},b_array, 'uniformoutput', false);
L_array=cellfun(@sys_prod,{TFopen},b_array, 'uniformoutput', false);
S_array=cellfun(@feedback,{tf(1,1)},L_array, 'uniformoutput', false);
[Gm,Pm,Wcg,Wcp] = cellfun(@margin,L_array);
disp(" ");
close all
figure 1
step(TFstage1,'r',A_array{:});
figure 2
bode(TFstage1,A_array{:})
figure 3
plot(C,Pm)
grid
xlabel('Compensation Capacitor, C (pF)');
ylabel('Phase Margin (deg)')
figure 4
step(A_array{C==50e-12},'r',A_array{C==100e-12},'b',A_array{C==200e-12},'g');
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 = A_array{C==100e-12};
bode(TFopen,'b',TFstage1,'g',TFcomp,'r');
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=20e-9;
R=1000;
TFsection=tf([1],[C*R,1])
TFfilter=series(series(TFsection,TFsection),TFsection);
disp("");
TFpreamp=series(series(TFcomp,TFfilter),TFcomp);
figure 1;
bode(TFpreamp,{1,1e5});
figure 2;
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;