Computational Complexity of Polyphase Resampling

Question

does someone know the approximate computational complexity of a polyphase resampling algorithm?

I am trying to find a good argument to prefer a polyphase resampling of Scipy over the FFT resampling method. In its documentation for resample_poly() it says:

This polyphase method will likely be faster than the Fourier method in scipy.signal.resample when the number of samples is large and prime, or when the number of samples is large and up and down share a large greatest common denominator.

My data meets those criteria. I resample 15 channels * 30k-50k samples from 10kHz to 300Hz and practice shows it is way faster than the FFT method.

Yet I would like to have a referable source on those criteria.

I found a few papers on the computational complexity of FFT algorithms to be around O(n*log(n)) with n=n_channel*n_samples.

Answer 1

Rational Resampling

10 kHz -> 300 Hz is a rational resampling with relatively benign factors:

$$ \frac{300\,\text{Hz}}{10\,\text{kHz}}=\frac{3\cdot10^2}{10^4}=\frac3{100}\text,$$

meaning that you'd interpolate by 3 and decimate by 100. That's exactly what you're achieving using resample_poly:

Polyphase Decimators

Remember that for decimation by 100 without aliasing, you need to filter your signal to $\frac1{100}$ of it's original bandwidth, which might take quite a lot of filter taps. Sadly, you throw away 99 of 100 samples that you filter in the following decimation step.

Using polyphase techniques, you can cleverly rearrange the $\frac1{100}$-band filter taps so that you don't even do these 99 calculations – and reduce the effort of your decimation filter by the decimation factor, here, 100.

This basically happens by taking your long filter and "deinterleaving" it to 100 disjunct sets, forming 100 polyphase components of your original filter. If you put in the input samples in all of these partial filters, and sum up the output samples of these component filters correctly, you get exactly the same filter. It's really just a way of "arranging" your filter differently. However, if you know which output samples you're going to throw away, you can omit pushing in each input sample to 99 of 100 polyphase filter components – so that you instead "round-robin" the input samples to each of the polyphase filter components. That saves 99% of calculations!

Polyphase resamplers

But what about the upsampling by 3? Well, another beautiful aspect of polyphase rational resamplers is that you can show that the order of "interpolate by 3 first, eliminate images by $\frac13$-band filter, then apply $\frac1{100}$-band filter to avoid aliasing, then throw away 99 of 100 samples" can be reordered to decimate-filter-interpolate. That's awesome because:

• the long filter runs at the lowest rate in the system
• you can omit the shorter filter, because it's far more relaxed than the longer one, anyway.

Why is that faster than the $\mathcal O(N+Nlog(N))$ of an FFT filter?

Aside from the aforementioned dramatic reduction in number of necessary operations, remember that even for the shorter polyphase component filters, you can apply FFT filters, if they are long enough so that this pays; in my experience, for rates as benign as yours, that's usually not the case. Compare a post by Tom Rondeau about comparison of FFT FIR filters to polyphase FIRs.

Generally, the beauty of FFT filters lies in the fact that you can reduce the number of operations you need to do to convolve a piece of signal with a filter tap vector. On the other hand, you buy that algorithmic elegance with totally non-linear memory accesses when computing that FFT. Now, if your filter is long, so must be your FFT, and your CPU can't keep all of both the signal vector and the tap vector in localized caches. Rule of thumb says that a single fetch from RAM is worth 200 to 500 Floating point operations^{[citation needed]} on modern x86 computers.

However, even if the amount of operations you need to handle is quite a lot larger in the non-FFT FIR filter, your memory access is (if your FIR implementation isn't totally braindead) inherently linear – and your CPU and its memory controller are really good at fetching linear memory, and even pre-fetching it. So, that works out in your favor quite a long time, until the complexity advantage of the FFT filter really beats the crap out of the naive FIR.

Another aspect of that linearity is easier vectorization using SIMD instructions – if you want to feel disappointed by scipy, install GNU Radio (easy if you're on current Linux distros or OS X, still doable on windows), run volk_profile (takes a while to find the optimum SIMD implementations for your machine), and then use the following flow graph to measure resampling throughput:

It does 5 million samples per second for me. Scipy's not getting quite close to that.