k=1000;
n=1e-9;
u=1e-6;

dp=2*pi;
f=100:10:5000;

%Crybaby circuit components
gf = -36  %forward gain of transistor stage
gp =  0   %gain of potentiometer

Cf = 10*n
Ci = 10*n
Rp = 33*k
Lp = 0.5
Ri = 68*k
Rs = 1.5*k

% Analog boiler-plate parameters
RpRi = Rp*Ri/(Rp + Ri);
f0 = 1/(2*pi*sqrt(Lp*Cf))
Q = RpRi*sqrt(Cf/Lp)

% DSP parameters
    M_PI = pi;
    fs = 44100
    w=dp.*f./fs;
    z1=e.^(-j.*w);
    z2=e.^(-2*j.*w);

    %Useful variables
    w0 = 2.0*M_PI*f0/fs;
    c = cos(w0);
    s = sin(w0);
    alpha = s/(2.*Q);
    
    %High Pass Filter coeffs
    b0h =  (1.0 + c)/2.0;
    b1h =   -(1.0 + c);
    b2h =  (1.0 + c)/2.0;
    a0h =   1.0 + alpha;
    a1h =  -2.0*c;
    a2h =   1.0 - alpha;
    
    %Band Pass Filter coeffs
    b0b =   Q*alpha;
    %b1b =   0.0;
    b2b =  -Q*alpha;
    a0b =   1.0 + alpha;
    a1b =  -2.0*c;
    a2b =   1.0 - alpha;
    
    %1-Pole HPF coeffs
    a1 = e^(-1/(Ri*Ci*fs));
    g = (1.0 + a1)*0.5;

%Transfer functions
    %1-Pole HPF
    Hin = g * (1 - z1)./(1 - a1.*z1);
    
    %Biquad BPF
    gbpf = 1/(Q*Q);
    Hbp = gbpf*(b0b  + b2b*z2)./(a0b + a1b*z1 + a2b*z2);

    %Biquad HPF
    Hhp = (b0h + b1h*z1 + b2h*z2)./(a0h + a1h*z1 + a2h*z2);
    
%semilogx(f, abs(Hbp) )

for gp = 0:0.1:1
    Hz = Hin.*( Rs/(Ri + Rs) + gf.*Hbp./(1 - gf.*gp.*Hhp) );
    semilogx(f, 20*log10(abs(Hz)) )
    %plot(f,abs(Hz))
    hold on
endfor

hold off


