2
$\begingroup$

I'm working on a Software Defined Radio project where I'd like to low-pass filter and decimate an analytical signal (IQ) sampled at 96ksps. Let's say the low-pass filter has a cutoff at 5kHz and I'd like to decimate by a factor 4 so that I have 24ksps out.

The idea is to perform the filtering using fast convolution using the overlap and save method as described in this article [pdf]:

http://www.3db-labs.com/01598092_MultibandFilterbank.pdf

I'm wondering if there are any pitfalls to my approach:

  1. Performing an N length FFT.
  2. Then doing an N length circular convolution (my multiplying with the FFT of my filter of length P).
  3. Then performing an N/4 IFFT back to decimate by 4 using the N/4 center taps of the forward FFT. Since my filter is a low-pass with a cutoff at 5kHz there should be very little energy outside the N/4 center taps of the FFT, and the P - 1 samples I need to discard should also lie outside the IFFT (if my filter is not too long).

EDIT:

This specific application is on a Raspberry Pi 3. After having given this some more thought I've realized it's not as clever as I first thought. I've stared too much on 0 centered FFTs and briefly forgot this is not the case. I would have to "remove" zeros in the middle of my FFT to make it shorter and then perform the IFFT.

What I will do is to do the fast convolution with an N length FFT and N length IFFT, and decimate when I copy samples to the output buffer.

$\endgroup$
8
  • 1
    $\begingroup$ hi! So, from an applicational perspective: Are you doing this on a microcontroller? Where does the requirement to optimize filtering a 96 kS/s convolution come from? And: if there is that requirement, why only decimate to 24 kS/s instead of, say, 6 kS/s? I'm currently trying to visualize the effect of your "middle fourth" IDFT; not quite sure this will work, but not convinced it doesn't, either. $\endgroup$ Jan 24, 2018 at 18:12
  • $\begingroup$ Edited the question to clarify. It's not as clever as I first thought. $\endgroup$ Jan 24, 2018 at 18:55
  • $\begingroup$ Ah, OK, is this is running on a fully-fledged OS (like Linux), you can just use existing signal processing software, e.g. GNU Radio, which comes with a decimating FIR Filter in FFT implementation. Again, there's really not much need to be overly efficient if this all happens at 96 kS/s: That's practically "nothing" computationally, and as long as your filter is not tens of thousands of taps long, your Pi (assuming a bit of NEON optimization, which GNU Radio has for most filtering applications) will easily manage, even with a direct decimating FIR (instead of fast convolution). $\endgroup$ Jan 24, 2018 at 20:44
  • 1
    $\begingroup$ @AlbinStigo: This is a well-established technique. See this article for a description of how it's done. $\endgroup$
    – Jason R
    Jan 24, 2018 at 21:28
  • 1
    $\begingroup$ @AlbinStigo: Sorry, I should have read the question more closely! $\endgroup$
    – Jason R
    Jan 25, 2018 at 20:25

1 Answer 1

6
$\begingroup$

I do a lot of decimation in the frequency domain. Little details are important.

I assume you already know the basic rules for fast convolution: the FFT length N is equal to the data blocksize L plus the length of the filter impulse response M minus 1. Each operation uses L samples of new data plus M-1 samples of data from the old block.

Ensure that the impulse response of your lowpass filter is shifted to the front of your time domain buffer AND properly windowed to M samples before you take the forward FFT to get the frequency domain representation of your filter. This keeps the result from wrapping around in the time domain when you take the inverse FFT. (Remember you're actually doing circular convolution when you want linear convolution.)

Kaiser is by far my favorite window because of its tuning knob. Use a large enough value to push the sidelobes down, as they will alias into your output. I typically construct a zero-phase brick wall in the frequency domain, take the inverse transform, apply the window and shift to the front of the buffer, then take the forward FFT.

Also make sure that the main lobe of your frequency response is still essentially zero past the Nyquist limit of your decimated output. I.e., don't make the Kaiser parameter too big.

If you're also doing frequency shifting by rotating the frequency bins, remember that you have to shift by a number of bins that corresponds to one data block size. For example, if L=M-1, then you can translate by any even number of bins.

Don't worry about FFT blocksizes that aren't powers of two. Just pick a convenient size that doesn't have any large prime factors and FFTW3 will perform well. I use L=M-1=3840 to decimate 192 kHz by 4:1 to 48 kHz with a block time of 20 ms. That's a FFT blocksize of L+M-1 = 7680 = 2^9 * 3 * 5.

$\endgroup$
3
  • $\begingroup$ Thanks for the comprehensive answer Phil! You can read more about our SDR at tujasdr.com if you are interested. 73 de SM6WJM $\endgroup$ Jan 20, 2019 at 17:23
  • $\begingroup$ Phil Karn ' May I ask you pls about decimation in frequency domain. And the concept of frequency shift by bins rotation ' how you can do it and why where I am so beginner. Could you pls put example for that Thanks alot $\endgroup$ Oct 17, 2019 at 8:40
  • $\begingroup$ I have asked a question about your answer [dsp.stackexchange.com/questions/88074/…. I would appreciate it if you could answer it. $\endgroup$
    – alireza
    Jun 5 at 14:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.