# Compare GSt's notch filter with LS's one in a frequency response plot.
# filename = path and file name of the respective IR sample .dat file
# midi_cutoff = GSt's abstract filter cutoff value (range of 0 to 127)
function res = plotall(filename, midi_cutoff)
# first calculate the real world center frequency (courtesy of Andreas Persson)
#c0 = midi_cutoff / 127;
#if (c0 < 0.5)
# c0 = c0 * 4826 - 1;
#else
# c0 = c0 * 5715 - 449;
#endif
#
# hmmm... or better least squares magic?
c0 = 39.08845189 * midi_cutoff - 0.03877717317 * midi_cutoff^2;
# now the approximation of GSt's notch filter ...
Bw = c0 * 1.357960934; # the so called "3dB bandwith" of our filter (the coefficient was calculated by "least squares matching" of all IR samples)
b = cos(c0 / 22050 * pi);
a = 1 + sqrt(1-cos(Bw/22050*pi)^2)/cos(Bw/22050*pi);
p = (1+a)/2;
num = [ p, -p*2*b, p ]; # numerator of our transfer function
den = [ 1, -b*(1+a), a ]; # denumerator of our transfer function
[fresp, frequencies] = freqz(num,den,44100);
frequencies = frequencies .* (22050 / pi);
title_redesign = [ ";c0=", num2str(c0), "Hz, Bw=", num2str(Bw), "Hz;" ];
# and here the respective IR sample for comparison
file = fopen(filename, "rb");
res = fscanf(file, "%f");
fclose(file);
y = fft(res);
len = floor((length(y) + 1) / 2);
y = y(2:len);
x = [1:len-1] / len * 22050;
title_file = [ ";", filename, ";" ];
# finally plot and compare our designed filter against the original one
axis([20, 22050, -90, 20]);
xlabel("Frequency [Hz]");
ylabel("Attenuation [dB]");
semilogx(frequencies, 20 * log(abs(fresp)), title_redesign, x, 28.7 + 20 * log(abs(y)), title_file);
endfunction;