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