# 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;
