Downsample a signal by a non-integer factor

Question

I have a signal with a sample rate of 8.9286 MHz and I want to downsample it to 500 kHz.

Since 8.9286 is not an integer multiple of 0.5 I can't simply decimate. Which downsample techniques are recommended to be used?

If I use this solution then I there will be a step to decimate by a factor of 647 $$(89286=2\times3\times23\times647),$$ which is prime, and would alias my signal in frequencies that are useful, so I can't filter them.

I could apply DFT to the original signal, then reconstruct in the time domain with the sample rate of 500 kHz, but that would take too much computational effort.

Is there a faster or simpler way to downsample it?

Answer 1

You need a resampler. There's different resamplers!

I'm not sure you're correctly interpreting the answer you cite: you don't just decimate by 647; you'd first (at least mathematically) upsample by some other factor. I'm not even quite where you got the 647 from, but an exact rational resampler would have to interpolate by 2500 and decimate by 44643. That's indeed a very unattractive resampler. 2500 and 44643 are coprime (or relatively prime), not prime in themselves; this is only important because if they weren't, you could cancel out factors.

Generally, you'd avoid resamplers where the interpolation or decimation sizes are so extreme, because you'll need to designa a $\frac1{44643}$-band filter, and that will lead to a very long filter.

Instead, you first try to bring the input rate "closer" to the output rate, e.g. by decimation-by-4, then apply an arbitrary resampler. Such arbitrary resamplers approximate the signal shape between samples at more or less arbitrary points. There's again many arbitrary resampler architectures, but the one I most commonly use is the polyphase filter bank arbitrary resampler, as e.g. implemented by the GNU Radio PFB arbitrary resampler block.

However, it's worth noting that in you fixed-input/output-rate scenario, it might be worth considering whether you actually need your output rate to be exactly 500 kHz. For example, assuming you use this for wireless communications with packets with a preamble. Then, using an upsampling by 20 and a downsampling by 357 will give you a rate that's only 4·10⁻⁴ larger than 500 kHz. If your packet is, say, only 256 samples long after the preamble, when you properly aligned timing in the preamble, the timing of your last payload sample will be 10% too early. You could even incorporate that knowledge, and retard your timing synchronization at the preamble by 5% of a symbol period, so that the maximum absolute timing error over the whole packet is but 5%, if 10% timing error increases your error vector magnitude too much. 20 up, 357 down is still not an easy rational resampling, but on my slightly older CPU core it runs, at an output rate of 5.1 MHz/s around 20% faster than the exact arbitrary resample, and about twice as fast as the exact rational resampler (2500/44643).

Notice that

1. even the exact arbitrary resampler is about five times as fast running on a single x86_64 CPU core than it needs to be to achieve 500 kHz, so if that's your platform, being "smarter" than just using an existing rational resampler implementation doesn't help (aside from maybe saving energy).
2. if you plan to implement this in hardware, check whether it would be faster / easier to hold a whole set of filter coefficient vectors in memory and evaluate a different filter for every output sample, or to have one very long filter. Use the polyphase arbitrary, or the rational resampling approach, accordingly.

Answer 2

If python is OK, you can just resample using an appropriate (approximate) interpolator:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.signal import butter, lfilter
T = 10000
t = np.linspace(0,1,T, endpoint=True)
x = np.sin(2np.pit*10)
ratio = 8.9286/0.5
t2 = np.linspace(0,1,int(T/ratio))
b, a = butter(6, 2*np.pi/ratio, btype='low', analog=False)
x_lpf = lfilter(b,a,x)
interpolated = interp1d(t, x_lpf, kind='cubic')
plt.plot(t,x)
plt.plot(t2,interpolated(t2),'.')