A question that pops up for many DSP-ers working with IIR and FIR filters, I think, is how to look at a filter\u2019s frequency and phase response. For many, maybe they\u2019ve calculated filter coefficients with something like the biquad calculator on this site, or maybe they\u2019ve used a MATLAB, Octave, Python (with the scipy library) and functions like freqz to compute and plot responses. But what if you want to code your own, perhaps to plot within a plugin written in c++?<\/p>\n
You can find methods of calculating biquads, for instance, but here we\u2019ll discuss a general solution. Fortunately, the general solution is easier to understand than starting with an equation that may have been optimized for a specific task, such as plotting biquad response.<\/p>\n
One way we could approach it is to plot the impulse response of the filter. That works for any linear, time-invariant process, and a fixed filter qualifies. One problem is that we don\u2019t know how long the impulse response might be, for an arbitrary filter. IIR (Infinite Impulse Response) filters can have a very long impulse response, as the name implies. We can feed a 1.0 sample followed by 0.0 samples to obtain the impulse response of the filter. While we don\u2019t know how long it will be, we\u00a0could take a long impulse response, perhaps windowing it, use an FFT to convert it to the frequency domain, and get a pretty good picture. But it\u2019s not perfect.<\/p>\n
For an FIR (Finite Impulse Response) filter, though, the results are precise. And the impulse response is equal to the coefficients themselves. So:<\/p>\n
For the FIR, we simply run the coefficients through an FFT, and take the absolute value of the complex result to get the magnitude response.<\/em><\/p>\n (The FFT requires a power-of-2 length, so we\u2019d need to append zeros to fill, or use a DFT. But we probably want to append zeros anyway, to get more frequency points out for our graph.)<\/p>\n Let\u2019s look for a more precise way to plot an arbitrary filter\u2019s response, which might be IIR. Fortunately, if we have the filter coefficients, we have everything we need, because we have the filter\u2019s transfer function, from which we can calculate a response for any frequency.<\/p>\n The transfer function of an IIR filter is given by<\/p>\n z0<\/sup><\/em>\u00a0is 1, of course, as is any value raised to the power of 0. And for normalized biquads, b0<\/sub><\/em> is always 1, but I\u2019ll leave it here for generality\u2014you\u2019ll see why soon.<\/p>\n To translate that to an analog response, we substitute ej\u03c9<\/sup><\/em> for z<\/em>, where \u03c9<\/em> is 2\u03c0*freq, with freq being the normalized frequency, or frequency\/samplerate:<\/p>\n Again, e0j\u03c9<\/sup><\/em> is simply 1.0, but left so you can see the pattern. Here it is restated using summations of an arbitrary number of poles and zeros:<\/p>\n For any angular frequency, \u03c9<\/em>, we can solve H(ej\u03c9<\/sup>)<\/em>. A normalized frequency of 0.5 is half the sample rate, so we probably want to step it from 0 to 0.5\u2014\u03c9<\/em> from 0 to \u03c0\u2014for however many points we want to evaluate and plot.<\/p>\n From that last equation, we can see that a single FOR loop will handle the top or the bottom coefficient sets. Here, we\u2019ll code that into a function that can evaluate either zeros (a<\/em> terms) or poles (b<\/em> terms). We\u001d\u2019ll refer to this as our direct evaluation function, since it evaluates the coefficients directly (as opposed to evaluating an impulse response).<\/p>\n You\u2019ve probably noticed the j<\/em>, meaning an imaginary part of a complex number\u2014the output will be complex. That\u2019s OK, the output of an FFT is complex too, and we know how to get magnitude and phase from it already.<\/p>\n Some languages support complex arithmetic, and have no problem evaluating \u201c For languages that don\u2019t, we can use Euler\u2019s formula, ejx<\/sup> = cos(x) + j * sin(x)<\/em>; that is, the real part is the cosine of the argument, and the imaginary part is the sine of it.<\/p>\n (Remember, j<\/em> is the same as i<\/em>\u2014electrical engineers already used i<\/em> to symbolize current, so they diverged from physicists and used j<\/em>. Computer programming often use j<\/em>, maybe because i<\/em> is a commonly used index variable.)<\/p>\n So, we create our function, run it on the numerator coefficients for a given frequency, run it again on the denominator coefficients, and divide the two. The result will be complex\u2014taking the absolute value gives us the magnitude response at that frequency.<\/p>\n Since we already had a precise method of looking at FIR response via the FFT\/DFT, let\u2019s compare the two methods to see how similar they are.<\/p>\n To use our new method for the case of an FIR, we note that the denominator is simply 1, so there is no denominator to evaluate, no need for division. So:<\/p>\n For the FIR, we simply run the coefficients through our evaluation function, and take the absolute value of the complex result to get the magnitude response.<\/em><\/p>\n Does that sound familiar? It\u2019s the same process we outlined using the FFT.<\/p>\n OK, we just showed that our new evaluation function and the FFT are equivalent. (There is a difference\u2014our evaluation function can check the response at an arbitrary frequency, whereas the FFT frequency spacing is defined by the FFT size, but we\u2019ll set that aside for the moment. For a given frequency, the two produce identical results.)<\/p>\n Now, if the direct evaluation function and the FFT give the same results, for the same frequency point, and the numerator and denominator are evaluated by the same function, by extension we could also get a precise evaluation by substituting an FFT process for both the numerator and denominator, and dividing the two as before. Note that we\u2019re no longer talking about the FFT of the impulse response, but the coefficients themselves. That means we no longer have the problem of getting the response of an impulse that can ring out for an unknown time\u2014we have a known number of coefficients to run through the FFT.<\/p>\n In general, the answer is our direct evaluation method. Why? We can decide exactly where we want to evaluate each point. That means that we can just as easily plot with log frequency as we can linear.<\/p>\n But, there may be times that the FFT is more suitable\u2014it is extremely efficient for power-of-2 lengths. (And don\u2019t forget that we can use a real FFT\u2014the upper half of the general FFT results would mirror the lower half and not be needed.)<\/p>\n We probably want to evaluate \u03c9 from 0 to \u03c0, corresponding to a range of half the sample rate. So, we\u2019d call the evaluation function with the numerator coefficients and with the denominator coefficients, for every \u03c9 that we want to know (spacing can be linear or log), and divide the two. For frequency response, we\u2019d take the absolute value (equivalently, the square root of the sum of the squared real and imaginary parts) of each complex result to obtain magnitude, and arc tangent of the imaginary part divided by the real part (specifically, we use the atan2 function, which takes into account quadrants). Note that this is the same conversion we use for FFT results, as you can see in my article, A gentle introduction to the FFT<\/a>.<\/p>\n For now, I\u2019ll leave you with some Python code, as it\u2019s cleaner and leaner than a C or C++ implementation. It will make it easier to transfer to any language you might want (Python can be quite compact and elegant\u2014I\u2019m going for easy to understand and translate with this code). Here\u2019s the direct evaluation routine corresponding to the summation part of the equation (you\u2019ll also need to \u201cimport numpy\u201d to have e<\/em> available\u2014also available in the math library, but we\u2019ll use numpy later, so we\u2019ll stick with numpy alone):<\/p>\n Again, we call this with the coefficients for each frequency of interest. Once for the numerator coefficients (the a<\/em> coefficients on this website, corresponding to zeros), once for the denominator coefficients (b<\/em>, for the poles\u2014and don\u2019t forget that if there is no b0<\/em>, the case for a normalized filter, insert a 1.0 in its place). Divide the first result by the second. Use use abs (or equivalent) for magnitude and atan2 for phase on the result. Repeat for every frequency of interest.<\/p>\n Here\u2019s a python function that evaluates numerator and denominator coefficients at an arbitrary number of points from 0 to \u03c0 radians, with equal spacing, returning arrays of magnitude (in dB) and phase (in radian, between +\/- \u03c0):<\/p>\n Here\u2019s an example of evaluating biquad coefficients at 64 evenly spaced frequencies from 0 Hz to half the sample rate (these coefficients are right out of the biquad calculator on this website\u2014don\u2019t forget to include b0 = 1.0):<\/p>\n Next up, a javascript widget to plot magnitude and phase of arbitrary filter coefficients.<\/p>\n The direct evaluation function performs a Fourier analysis at a frequency of interest. For better understanding, reconcile it with the discrete Fourier transform described in A gentle introduction to the FFT<\/a>. In that article, I describe probing the signal with cosine and sine waves to obtain the response at a given frequency. Look again at Euler\u2019s formula, which shows that ej\u03c9<\/sup><\/em> is cosine (real part) and sine (imaginary part), which the article alludes to this under the section \u201cGetting complex\u201d. You should understand that the direct evaluation function presented here could be used to produce a DFT (given complete evaluation of the signals at appropriately spaced frequencies). The main difference is that for this analysis, we need not do a complete and reversible transform\u2014we need only analyze frequency response values that we want to graph.<\/p>\n","protected":false},"excerpt":{"rendered":" A question that pops up for many DSP-ers working with IIR and FIR filters, I think, is how to look at a filter\u2019s frequency and phase response. For many, maybe they\u2019ve calculated filter coefficients with something like the biquad calculator … Continue reading Plotting the filter precisely<\/h3>\n
\\(H(z)=\\frac{a_{0}z^{0}+a_{1}z^{-1}+a_{2}z^{-2}…}{b_{0}z^{0}+b_{1}z^{-1}+b_{2}z^{-2}…}\\)<\/h3>\n
\\(H(e^{j\\omega})=\\frac{a_{0}e^{0j\\omega}+a_{1}e^{-1j\\omega}+a_{2}e^{-2j\\omega}…}{b_{0}e^{0j\\omega}+b_{1}e^{-1j\\omega}+b_{2}e^{-2j\\omega}…}\\)<\/h3>\n
\\(H(e^{j\\omega})=\\frac{\\sum_{n=0}^{N}a_{n}e^{-nj\\omega}}{\\sum_{m=0}^{M}b_{m}e^{-mj\\omega}}\\)<\/h3>\n
Coding it<\/h3>\n
e**(-2*j*0.5)<\/code>\u201d\u2014either directly, or with an \u201cexp\u201d (exponential) function. It\u2019s pretty easy in Python, for instance. (Something like, coef[idx] * math.e**(-idx * w * 1j)<\/code>, as the variable idx<\/code> steps through the coefficients array.)<\/p>\nRevisiting the FIR<\/h3>\n
And back to IIR<\/h3>\n
Which is better?<\/h3>\n
An implementation<\/h3>\n
\\(magnitude:=\\left |H \\right |=abs(H)=\\sqrt{H.real^2+H.imag^2}\\)<\/h5>\n
\\(phase := atan2(H.imag,H.real)\\)<\/h5>\n
import numpy as np\r\n\r\n# direct evaluation of coefficients at a given angular frequency\r\ndef coefsEval(coefs, w):\r\n res = 0\r\n idx = 0\r\n for x in coefs:\r\n res += x * np.e**(-idx * 1j * w)\r\n idx += 1\r\n return res<\/pre>\n
# filter response, evaluated at numPoints from 0-pi, inclusive\r\ndef filterEval(zeros, poles, numPoints):\r\n magdB = np.empty(0)\r\n phase = np.empty(0)\r\n for jdx in range(0, numPoints):\r\n w = jdx * math.pi \/ (numPoints - 1)\r\n resZeros = coefsEval(zeros, w)\r\n resPoles = coefsEval(poles, w)\r\n\r\n # output magnitude in dB, phase in radians\r\n Hw = resZeros \/ resPoles\r\n mag = abs(Hw)\r\n if mag == 0:\r\n mag = 0.0000000001 \u00a0# limit to -200 dB for log\r\n magdB = np.append(magdB, 20 * np.log10(mag))\r\n phase = np.append(phase, math.atan2(Hw.imag, Hw.real))\r\n return (magdB, phase)<\/pre>\n
zeros = [ 0.2513643668578741, 0.5027287337157482, 0.2513643668578741 ]\r\npoles = [ 1.0, -0.17123074520885395, 0.1766882126403502 ]\r\n\r\n(magdB, phase) = filterEval(zeros, poles, 64)\r\n\r\nprint(\"\\nMagnitude:\\n\")\r\nfor x in magdB:\r\n print(x)\r\n\r\nprint(\"\\nPhase:\\n\")\r\nfor x in phase:\r\n print(x)<\/pre>\n
Extra credit<\/h3>\n