% Calculate complex response of a two-pole low pass filter.
function response = twopole(Q, freq, fc, fn)
% Q = quality factor
% freq = array of frequencies
% fc = critical frequency (optional, default=1)
% fn = f0/fc (optional, calculated to put -3dB at fc if not supplied)

if not(exist('fc')) % Implement default critical frequency
   fc = 1.0;
endif

omega = 2*pi*freq; % angular frequency
omega_c = 2*pi*fc;

if not(exist('fn')) % Find relative location of the -3dB point
   b = 1 - 1/(2*Q^2);
   fn = 1/sqrt(b + sqrt(b^2 + 1));
endif

omega0 = omega_c * fn;

response = 1./(1 + j.*omega./Q./omega0 - (omega./omega0).^2 );
