Fast linear-phase filtering

[Rolfe]

I need to band-pass filter a large quantity of data in real time without introducing frequency-dependent group delay. I've previously tried generating FIR filter coefficients using the Parks-McClellan VI and running that along with my data (in chunks of around 32000 samples) through the convolution VI (with the default frequency-domain setting), but I found this to be too slow. The options I'm considering at this point are

a) writing my own convolution VI which uses an overlap-add method (ala Oppenheim and Schafer, pg 558) so that the FFT computations are effectively linear time instead of O(n log n). (Specifically, filter small sections of the input signal by applying a fourier transform, multiplying with the transform of the FIR filter, inverse transforming (being careful to avoid time-aliasing). Then recombine the sections by adding the FIR tail from one filtered section into the signal for the next section before filtering)

b) running my signal through an IIR filter (eg Chebyshev), time-reversing the output, running the time-reversed filtered signal through the same filter again, time-reversing the output of the second filter. The second filtering operation reverses all group delays from the first filter so the total operation is zero-phase. That is, I would like to copy the matlab function filtfilt, but the labview documentation suggests that some tricks are required to handle the edges of the signal

c) try to construct an all-pass filter with group delay which is the additive inverse (plus a constant) of the group delay of my original IIR band-pass filter.

While I would love any general advice, I was wondering in particular about the implementation details of the LabVIEW convolution vi. Specifically, does it already use an overlap-add method when operating in frequency-domain mode? Alternatively, have any of you written a VI which will double-filter using an IIR filter to eliminate the non-linear phase? Perhaps someone knows a good reference for an algorithm for constructing all-pass filters with specified phase characteristics?

Thanks for your help,
Jason

[LocalDSP]

Jason,

These are very interesting questions and we will come back to you shortly with more answers. In the mean time let me give you some quick info:

The actual LabVIEW Convolution is, as you probably noticed, a 'one-shot' function that will not work properly in a filtering of consecutive data blocks situation since it will return a data set longer than your input record that includes 'settling' and 'un-settling' of your filter. To get it work correclty is not that difficult and we do have example code that I can send to you.

The convolution in Frequency-mode is based on simple FFT/IFFT operations on the input signals that are appropriately zero-padded to avoid the effects of circular convolution.

Running an IIR filter forward and backward is a
very good and efficient way to achieve linear phase but the method conflicts with your real-time requirements since you can't run the 'backward' filter until you are completely done with your acquisition (due to the Infinite Impulse Response nature of the filter). I guess it should be possible to achieve a reasonable result by cutting the length of your impulse response and work on a couple of acquisition blocks simultaneously, but the management of the different signal pieces becomes a complicated story (but not impossible)

Another technique you may consider is a combination of IIR and FIR filters, for example using a Bessel IIR type IIR filter with good phase linearity extended with a low-order FIR filter to correct the passband magnitude characteristic (PMC can design the needed FIR filter for you). There are different ways to do that and you can find ideas in the literature.

More to come...
regards, NI-DSP/Analysis

[LocalDSP]

PS: What version of LabVIEW are you using?

[LocalDSP]

Jason,

Here are a couple of options and tools that may work for you. In the attached VI library I have included two examples that:

1 - Perform a continuous convolution (frequency domain mode) with the Impulse Response (IR) of an FIR filter designed with Park-McClellan. This VI is optimized and should run faster than before, and since it also works correctly on consecutive data blocks it may solve your problem.

2 - Perform a two steps filtering by first using an IIR filter and then performing a convolution with a truncated and reversed version of the IR of the IIR filter. Depending on the needed truncation length this may increase performances compared to a direct FIR.

Let me know if this work for you otherwise we'll have to look at ot
her methods.

[Rolfe]

Thanks for the response. I'm using LabVIEW 7.1, so I'm pretty sure that the default behavior for convolution.vi is to do frequency-domain convolution. Is that incorrect? Is the frequency-domain convolution VI you just sent me more efficient than the default version?

I'm already performing the overlap-and-add step which is the basis of the continuous-convolution VI included in your last message, so I'm basically already using your first suggestion. I imagine I could improve the efficiency slightly if I performed the Fourier transform on the kernel only once rather than each time a chunk of data is filtered, but I don't think that would yield substantial savings, since FFT is O(n log n) in the size of the input, and the kernel I'm using (about 1000 taps) is an order of magnitude shorter than the data signal I'm using(about 16000 samples for each application of the filter).

The real issue is that IIR filtering seems to run about 5-10 times faster than FIR filtering. While your second suggestion would certainly work, my intuition is that it would be little faster than your first suggestion since it depends on FIR filtering for the second step. Notably (but not surprisingly), the speed of frequency-domain convolution is not substantially affected by changing the number of taps in the kernel. I've included a simple little VI which does both FIR (Parks-McClellan) and IIR (Chebyshev) filter. By profiling the operation of this VI, you can see the disparity in running time between FIR and IIR filtering.

At the point, double-IIR filtering seems like the best option. I haven't tested it yet, but I would think that the transients would moslty die out after a number of samples beyond the end of a data chunk comparable to the order of the filter. Given that I only need a 7th-order Chebyshev filter to get a good frequency response curve, I think I should be able to do some black magic to smooth over the transitions between data blocks in real time.

However, I'm concerned that actually time-reversing a signal after filtering will require me to copy the signal from one location to another, and thus will take a long time in LabVIEW. Does Reverse 1D Array work in some clever way that doesn't require actually moving the data? To avoid this problem, I figure I can just reimplement IIR Cascade Filter.vi and modify it to read the array backwards rather than forwards. Does this sound reasonable?

[LocalDSP]

The Convolution VI default behavior is 'Frequency Domain' but the VI has an input selector that you can change to "direct" that is classic time domain convolution. The Frequency Domain Convolution VI I posted yesterday was faster than the one shipping with 7.0, but the performances were significantly improved for 7.1 so please discard that VI.

You are right in the fact that a shorter kernel will not speed up the frq.domain convolution since the FFT size is mainly decided by the long input signal. However when the kernel gets short enough you hit the trade-off point where time domain convolution gets faster than frq. domain. I have attached a benchmark VI where you can select the kernel size and the convolution mode. With a
data set of 32.000 samples the trade off point on my machine corresponds to a kernel size of 250 coefficients. For your one shot FIR filtering with 1000 coefficient the frq. domain convolution is therefore the right choice but in the combined IIR+FIR example I posted yesterday, the correction FIR size was only 20 coefficients and, in this case, the time domain convolution is approximatelly 5 times faster than frq. domain.

But all together using your filter specifications and implementing the IIR+FIR (time domain convolution) you can get a factor of two in speed improvement compared to frq. domain. I have also attached a revised version of that VI adapted to your specifications.

Finally to answer your Reverse 1D Array question. It's a little tricky, because internally the reverse array function will not allocate a new buffer (just reverse the pointer values), but as soon as you connect the output to a sub-VI or a Call Library Function the data will be copied. So you can unfortunat
ely not get around this.