EPx professional blog and repository for braindumps

2010/01/18

IIR filter, in Python (as usual)

And now, for completeness, let's have an IIR (infinite Impulse Response) filter. Those digital filters calculate the output signal based not only on past input signals, but also on past OUTPUT signals. In other words, there is a feedback component on IIR output. This filter may have output forever, even if the input signal is cut, hence the "infinite" quality.

An IIR filter can be expressed by a difference equation, which defines output in terms of previous outputs, as well as inputs. A very popular example of difference equation is the Fibonacci sequence:


Fib(n) = Fib(n-1) + Fib(n-2)


In digital filtering, we are more interested in difference equations that are stable, that is, tend to zero as the input fades away. Since a functional IIR filter can be defined by a difference equation of just two terms, like the Fibonacci series above, it can do much more using much less math. In this sense, IIR is much more powerful than FIR.

The following example is a kind of resonator that enhances frequency response around 2000Hz (any other value could be chosen just by changing the source code). The "r" factor is the resonance power. If you make r=1, it becomes an oscilator tuned to 2000Hz. If you make r > 1, it becomes exponential and useless.


#!/usr/bin/env python

import wave, struct, math
SAMPLE_RATE = 44100

original = wave.open("NubiaCantaDalva.wav", "r")
filtered = wave.open("NubiaFilteredIIR.wav", "w")
filtered.setnchannels(1)
filtered.setsampwidth(2)
filtered.setframerate(SAMPLE_RATE)

# IIR filter coefficients
freq = 2000 # Hz
r = 0.98
a1 = -2.0 * r * math.cos(freq / (SAMPLE_RATE / 2.0) * math.pi)
a2 = r * r
filter = [a1, a2]
print filter

n = original.getnframes()
original = struct.unpack('%dh' % n, original.readframes(n))
original = [s / 2.0**15 for s in original]

result = [ 0 for i in range(0, len(filter)) ]
biggest = 1
for sample in original:
for cpos in range(0, len(filter)):
sample -= result[len(result) - 1 - cpos] * filter[cpos]
result.append(sample)
biggest = max(abs(sample), biggest)

result = [ sample / biggest for sample in result ]
result = [ int(sample * (2.0**15 - 1)) for sample in result ]
filtered.writeframes(struct.pack('%dh' % len(result), *result))


You can see that our IIR filter can be directly derived from frequency and resonance factor, and it has only two weights (a1 and a2), so it needs only 2 multiplications per sample, while in FIR we needed dozens, or hundreds.

Another thing you may notice is the "biggest" variable. IIR filters are quite the opposite of FIR filters; they can only increase frequency response in a certain range, while FIR can only decrease it. So, in IIR filter, we let it increase (a lot) the response around 2000Hz (above the clipping limits) and then scale all values back to avoid clipping. The net effect is of course damp all frequencies except those around 2000Hz.

The IIR filter coefficients are NOT the impulse response like it were in FIR; so you can NOT find the filter weights using an inverse FFT as we did in FIR. The impulse response of an IIR filter can be found by FFT, but this is not directly useful for actual time-domain processing.

In fact, IIR filters are VERY difficult to calculate and prone to instability problems due to feedback. Since at least the two-order IIR filter is easy to make, as we did above, there are methods to make higher-order IIR filters as a bunch of two-order IIR sub-filters, connected in series or parallel.

Another way (that I have heard of) to design IIR filters is by mimicking some well-known analog filter design. In fact, all analog filters have "infinite response" and rely on feedback (even though analog circuits have continous output), so they match IIR features well.

All software-defined filters I have seen are based on FIR. Maybe embedded things like guitar pedals are using IIR, but I can't tell. Since I don't possess the needed math to realize more complex IIR filters, and the FFT->FIR weights method worked so well, I'd go FIR if I needed to design a filter right now...

0 comentários:

Postar um comentário