Inverting amplifier filter analysis

filter_test.m

% Here's an attempt to sketch the filter toplogy:
%
% in      R1A       C1            R2A       C2
% >----+--/\/\/\----||---+--+--+--/\/\/\----||---+---+----> out
%      |                 |  |  |                 |   |
%      +------/\/\/\-----+  |  +------/\/\/\-----+   |
%             R1B           |         R2B            |
%                           +--|-\                   |
%                              |  \__________________|
%             R_BIAS           |  /
%   GND|------/\/\/\-----------|+/ OPAMP
%
% The transfer function of this topology is:
%
%         R2B   w_p1 w_p2   (s - w_z1) (s - w_z2)
% H(s) = ----- ----------- -----------------------
%         R1B   w_z1 w_z2   (s - w_p1) (s - w_p2)
%
% w_z1 = 1/([R1A + R1B] C1)     w_p1 = 1/(R1A C1)
% w_z2 = 1/(R2A C2)             w_z2 = 1/([R2A + R2B] C2)
%
% where w ("omega") are the angular frequencies of 
% the poles and zeros.  Note that the DC gain is -R2B/R1B.
%
% R_BIAS = R1B || R2B
%
% The purpose of R_BIAS is to balance the DC source impedance seen
% by each of the opamp inputs (see Horowitz & Hill, page 195).  It
% does not affect the transfer function, but does add noise.

f = logspace(log10(0.1), log10(5e6), 100);
omega = 2*pi*f;

parallel  = @(x,y) 1./(1./x + 1./y);
series    = @(x,y) x + y;
capacitor = @(c) 1./(j*omega*c);
resistor  = @(r) r * ones(size(omega));

% Physical constants
k_B = 1.3806503e-23; %  Boltzman's constant [m^2 kg / (s^2 K)]

% Operating conditions
T = 27 + 273.15;     %  Operating temperature [Kelvin]

% Opamp noise parameters (from LT1125 datasheet)
i_n = sqrt(2) * 1.3e-12; % current noise density [A/rtHz] at 10Hz
e_n = sqrt(2) * 3.0e-9;  % voltage noise density [V/rtHz] at 10Hz

%i_n = 1.7e-12;  % OP27
%e_n = 3.8e-9;   % OP27

R1A = 1000;
C1  = 4e-6;
R1B = 5000;

R2A =  1000;
C2  =  0.1e-6;
R2B =  5000;
C2B =  0;%3.3e-12;

R_BIAS = 2500; %parallel(R1B, R2B)

omega_z1 = 1/(C2 * R2A)
omega_p1 = 1/(C1 * R1A)
omega_z2 = 1/(C1 * (R1A + R1B))
omega_p2 = 1/(C2 * (R2A + R2B))

Z1 = parallel(resistor(R1B), series(resistor(R1A), capacitor(C1)));
Z2 = parallel(parallel(resistor(R2B), series(resistor(R2A), capacitor(C2))), ...
              capacitor(C2B));

G1 = -Z2./Z1;

% compute the noise
i_A1_sq = i_n^2 + 4*k_B*T ./ real(Z1) + 4*k_B*T ./ real(Z2);    % current noise
e_A1_sq = abs(1 + Z1./Z2).^2 .* e_n^2 + abs(Z1).^2 .* i_A1_sq; % voltage noise
e_A1_sq = e_A1_sq + abs(1 + Z1./Z2) * 4*k_B*T * R_BIAS;     % FIXME: Is this correct?

% Make a Bode plot
subplot(3,1,1);
semilogx(f, 20*log10(abs(G1)));
grid on;
title('Filter Module 1');
ylabel('gain [dB]');
axis tight;

% Display the poles and zeros
line([omega_z1 omega_z1]/(2*pi), get(gca, 'YLim'), 'color', 'black');
line([omega_z2 omega_z2]/(2*pi), get(gca, 'YLim'), 'color', 'black');
line([omega_p1 omega_p1]/(2*pi), get(gca, 'YLim'), 'color', 'red');
line([omega_p2 omega_p2]/(2*pi), get(gca, 'YLim'), 'color', 'red');

legend('response', ...
    sprintf('zero at %0.0f Hz', omega_z1/(2*pi)), ...
    sprintf('zero at %0.0f Hz', omega_z2/(2*pi)), ...
    sprintf('pole at %0.0f Hz', omega_p1/(2*pi)), ...
    sprintf('pole at %0.0f Hz', omega_p2/(2*pi)));

subplot(3,1,2);
semilogx(f, angle(-G1)*180/pi);
axis tight;
ylabel('phase [degrees]');
grid on;
set(gca, 'YLim', [-180 180], 'YTick', 45*(-4:4));


% Display the noise

subplot(3,1,3);
semilogx(f, sqrt(e_A1_sq));
hold all;
plot(f, sqrt(e_A1_sq).*abs(G1));
hold off;
legend('input-referred noise','output-referred noise');
grid on;
title('noise');
ylabel('volts / rtHz');
axis tight

[y, yi] = min(abs(f - 100));
G100 = abs(G1(yi));
fprintf('gain at 100Hz:   %f (%f dB)\n', G100, db(G100));
fprintf('inoise at 100Hz: %f nV/rtHz\n', sqrt(e_A1_sq(yi))*1e9);
fprintf('onoise at 100Hz: %f nV/rtHz\n', sqrt(e_A1_sq(yi))*1e9 * G100);