N = 1000;


frequency_range = 0:N - 1;
zfrequency_range = exp(i * pi * frequency_range / N);


lower_zcutoff = 0.01;

zorder = 2;
figure_number = 1;


[bfzeros, bfpoles, bfgain] = butter(zorder, lower_zcutoff);

znum = ones(1, length(frequency_range));
for zero_number = 1:length(bfzeros)
  znum = znum .* (bfzeros(zero_number) - zfrequency_range);
endfor

zdenom = ones(1, length(frequency_range));
for pole_number = 1:length(bfpoles)
  zdenom = zdenom .* (bfpoles(pole_number) - zfrequency_range);
endfor

frequency_response = bfgain * znum ./ zdenom;

figure(figure_number++);
semilogx(frequency_range(2:end), abs(frequency_response(2:end)));
title("magnitude response");

min_phase = -imag(hilbert(log(abs(frequency_response))));


figure(figure_number++);
semilogx(frequency_range(2:end), 180 * unwrap(arg(frequency_response(2:end)), pi) / pi);
title("measured phase response in degrees");

figure(figure_number++);
semilogx(frequency_range(2:end), 180 * unwrap(min_phase(2:end), pi) / pi);
title("min phase response in degrees");