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");


phase_response = unwrap(arg(frequency_response), pi);

min_phase = unwrap(-imag(hilbert(log(abs([frequency_response, fliplr(frequency_response(2:end))])))), pi); # thanks to Ben Abbot for 'fliplr' part


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

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

figure(figure_number++);
semilogx(frequency_range(2:end), 180 * (phase_response(2:end) - min_phase(2:numel(frequency_range))) / pi);
title("phase responses difference in degrees");