Draw a frequency response graph from a simple equation

Question

How is it possible to draw an audio frequency response graph (x-axis : freq from 20 to 20 000hz, y-axis : dB) from a simple equation like

y[n] = x[n] - x[n-486]

(Example: with this equation I get a nice comb filter that removes frequency 197,5Hz and all its harmonics, for an audio file with sampling rate 96 Khz)

How to turn such an equation into a graphical interpretation like this http://downloadfreesamples.com/wp-content/uploads/2013/08/EQ-Four-CIS-DSP_3.jpg ? Then I could see how deep / steep the filtering is.

Thank you

(PS : I continue to read references and online books about filter design that were given to me in previous questions, thanks for these references by the way.)

Answer 1

The tool of choice here is the z-transform, which is applied to your equation. It transforms sequences of complex numbers to functions in the complex plane. Its most useful properties are that it takes delays to powers of the argument z, and for z = exp(i omega) we get the fourier transform of the transformed sequence. Look up the details, I'll just show you how it works here.

y[n] = x[n] - x[n-486]  ---> Y(z) = X(z) - X(z)*z^(-486)

The equation in z-domain factors on the right hand side to give Y(z) = X(z) ( 1 - z^(-486) ).

Since Y(z) is the output and X(z) is the input, for linear systems we define the transfer function to be Y(z)/X(z) = H(z) and the above equation becomes

H(z) = Y(z)/X(z) = 1 - z^(-486)

Now if we substitute z = exp(i omega) we get the frequency response of your system, which is

H_f(omega) = 1 - exp(i omega)^(-486)

and with the properties of exp() it simplifies to

H_f(omega) = 1 - exp(- 486 i omega)

If you need actual frequency f instead angular frequency omega, just substitute omega = 2 pi f/fs, with sampling rate fs.

Finally, you'll probably want to know just the magnitude of the response, so look at abs(H_f). In this specific case it's very simple to calculate because abs(H_f) = sqrt( H_f * conj(H_f) ) and conj(H_f) = 1 - exp(+ 486 i omega).

As requested, some code. Here's what you can do in matlab/octave to get the magnitude frequency response with a logarithmic frequency axis and magnitude axis in dB.

% Create a logarithmically spaced frequency axis with range 20Hz to 20kHz

fAxis = logspace( log10( 20 ) , log10( 20000) , 10000);

% calculate the complex frequency response at these frequencies

% Set SampleRate before calling this, like SampleRate = 44100;

Hf = 1 - exp(1i*486*2*pi*fAxis/SampleRate);

% plot the magnitude spectrum

semilogx(fAxis,10*log10(Hf.*conj(Hf)));

Answer 2

Factor 1-x^486 in the complex domain. You will get a bunch of zeros on the unit circle. Calculate the product of the distances to all these zeros as you travel around the unit circle from pi*(20/96000) to pi*(20/96). Plot.

Answer 3

I finally found a solution (probably not the best, but ...)

1) Generate a 10 seconds white noise (which has flat frequency spectrum on 10 seconds average)

2) Apply the filter

3) Display the spectrum in a Sound editor software (e.g. Sound forge) with averaging on the 10 seconds.