How can I vectorize the computations for a first-order recursive filter?

I have a simple single pole low pass filter (for parameter smoothing) that can be explained by the following formula:

$$y[n] = (1-a) y[n-1] + a x[n]$$

The architecture that I'm using has access to single-instruction, multiple-data (SIMD) instructions that can perform multiple vectorized calculations in parallel. I would like to take advantage of this capability, but I'm not sure how to do that for a recursive filter like this. The problem is that every computation needs a previous result.

migrated from stackoverflow.comJan 5 '12 at 23:46

This question came from our site for professional and enthusiast programmers.

• Could someone please clarify why this was closed as "off topic" ? – Paul R Jan 6 '12 at 16:38
• The question overlaps between here and Stack Overflow. The original question asked specifically for how to implement it using ARM NEON extensions. The question will live on both sites; it has been edited here to make it more of a theoretical discussion on structuring the problem to take advantage of parallelism. – Jason R Jan 6 '12 at 16:41
• @PaulR I requested it to be migrated here yesterday, but phonon felt strongly that it belonged on SO, and not here. I'll admit, I don't know about ARM NEON in particular, and he's probably right and I respect his judgement. See this meta question. My deleted answer basically said that I agree with Jason R and that users should flag such questions for migration – Lorem Ipsum Jan 6 '12 at 16:42
• Thanks both for the clarification - I've been doing my best on SO to get DSP-related questions migrated here in order to raise our profile, so I was concerned as to whether this was the right thing to do when this question got closed - glad to see it's open again now in a more general form. – Paul R Jan 6 '12 at 17:19

Assuming that you perform vector operations $M$ elements at a time, you can unroll the difference equation by a factor of $M$ pretty easily for the simple single-pole filter. Assume that you have already calculated all outputs up to $y[n]$. Then, you can calculate the subsequent ones as follows: \begin{align} y[n+1] &= (1-a)y[n] + ax[n+1] \\ y[n+2] &= (1-a)y[n+1] + ax[n+2] \\ &= (1-a)((1-a)y[n] + ax[n+1]) + ax[n+2] \\ &= (1-a)^2y[n] + a(1-a)x[n+1] + ax[n+2] \\ &\vdots \end{align}

In general, you can write $y[n+k]$ as:

$$y[n+k] = (1-a)^ky[n] + \sum_{i=1}^k a(1-a)^{k-i}x[n+i]$$

For each sample index $n+k$, this looks like an FIR filter with $k+1$ taps: one tap multiplies the last filter output $y[n]$, and the other $k$ taps multiply the filter inputs $x[n+1], \ldots , x[n+k]$. The nice thing is that the coefficients used for all of these taps can be precalculated, allowing you to unroll the recursive filter into $M$ $M+1$-tap parallel non-recursive filters (these filters compute the output samples $y[n+1], \ldots , y[n+M]$), updating the feedback term every $M$ output samples. So, given an initial condition $y[n]$ (which is assumed to be the last output calculated on the previous vectorized iteration), you can calculate the next $M$ outputs in parallel.

There are some caveats to this approach:

• If $M$ becomes large, then you end up multiplying a bunch of numbers together in order to get the effective FIR coefficients for the unrolled filters. Depending upon your number format and the value of $a$, this could have numerical precision implications.

• Also, you don't get an $M$-fold speedup with this approach: you end up calculating $y[n+k]$ with what amounts to a $k$-tap FIR filter. Although you're calculating $M$ outputs in parallel, the fact that you have to do $k$ multiply-accumulate (MAC) operations instead of the simple first-order recursive implementation diminishes some of the benefit to vectorization. The non-vectorized approach uses 2 MACs per output, so you need $2M$ operations to calculate $M$ outputs. The vectorized scheme computes the $M$ outputs at once, requiring $M+1$ MACs in the process. So, the reduction in operations can be expressed as a function of $M$ as:

$$R = \frac{M+1}{2M} = \frac{1}{2}\left(1 + \frac{1}{M}\right)$$

So, even with $M$ very large, you can get a maximum of a 50% reduction in the amount of computations using this method. For the common values of $M=4$ and $M=8$, the reduction is 37.5% and 43.75%, respectively. However, in addition to pure reduction in the number of operations, you may gain additional performance benefits due to the different memory access pattern employed by the unrolled approach; for the simpler implementation, you may run into delays due to the dependency of each output sample $y[n]$ on the previous sample $y[n-1]$. This effect is obviously very platform-dependent.

In general, you can only vectorize completely independent sets of computations. But in your IIR low pass, every output is dependent on another (except the 1st), so vectorization is not possible.

If your variable "a" is large enough that (1-a)^n quickly decays to below your desired noise floor or allowed error, you could substitute a short FIR filter approximation for your IIR, and vectorize that convolution instead. But that's not likely to be faster.

You can only really vectorize this effectively if you have more than one signal to which you wish to apply the same filter, e.g. if it's a stereo audio signal then you can process the left and right channel in parallel. Four or eight channels in parallel would obviously be even better.

How about expanding equations to 4 steps and use matrix multiplication? a is constant so one matrix may be precalculated

It is most efficient to vectorise the filtering of several independent streams in parallel.

If you just have a single long stream you could also split the stream into several overlapping sections and filter these as if they were independent streams.

You want an overlap because the first few samples of each stream will be incorrect. The amount of samples you need to discard will depend on the value of a and the desired precision. You will want to discard about n samples where (1/a)(1-a)**n < eps in order that the result will be accurate to eps*max(|x[i]|).

For example, if a=0.1, then in 128 samples the error caused by (1/a)(1-a)^n would be less than the LSB in a 16bit integer, while for a=0.9 you would only need to discard about 6 samples.