% Calculate complex response of a multi-stage lowpass filter. Each stage
% is a 2-pole filter.
function response = multistage(freq, Qarr, FNarr, fc)
% freq = array of frequencies for which the response will be calculated.
% Qarr = array of quality factors (one for each stage)
% FNarr = f0/fc
% fc = critical frequency (optional, default=1)

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

Nstages = numel(Qarr);
if numel(FNarr) != Nstages
   error('Qarr and FNarr do not have the same number of elements, but they should.')
endif

response = 1; %initialize
for stage = 1:Nstages
   response = response .* twopole(Qarr(stage), freq, fc, FNarr(stage));
endfor