# Draw a frequency response graph from a simple equation

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.)

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

• Thanks! If I just want to draw the function with x=20..20000, what should I plot here ? plot(...,x=20..20000) (here with Maple, I will adapt from any other language) (I just need the frequency response graph, like in the screenshot I have put in my own answer in this topic here) – Basj Oct 19 '13 at 22:27
• Jazzmaniac, How to find the magnitude spectrum of some biquad filters in cascade ? – Basj Nov 9 '13 at 22:01

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.

• What language do you usually use for this? I'll try to code this with Scilab or R or something else... – Basj Oct 19 '13 at 22:11

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.

• Well it gives the "rough" idea, but it unfortunately doesn't work very precisely. I tried this method with a Notch (-52 dB attenuation on a Pure sine wave of the centre frequency!) and on the chart generated by this method we can only see 10 dB attenuation... So this method only gives a "rough idea", but is not precise! – Basj Oct 19 '13 at 20:55