How to do high-resolution FFT on just the lower frequencies in a signal?

Question

Description of the data and problem:

I have a signal sampled at 1000 Hz. I'm low-pass filtering it at 120 Hz, and want to make spectrograms of the frequencies below this threshold. I'm using the scipy functions fftpack.fftfreq to get the Fourier coefficients, then fftpack.fft to do the actual transform.

The signal is quite long, and I want the spectrograms to be about 5 seconds in length using 50 millisecond windows. I also filter out negative frequencies for plotting. Given the number of samples in this window (50), I get 50 Fourier coefficients. However, half of these are negative, and another portion goes between 120 Hz and 500 Hz (naturally, because the sampling rate is 1000 Hz).

This leaves me with pretty low-resolution spectrograms. I only have about 6 blocks in the frequency area of interest (0, 20, 40, 60, 100, and 120 Hz), then 18 blocks showing low activity in the frequencies up to 500 Hz, since they were filtered out.

Question:

How could I, for instance, do the FFT specifying something like 50 frequency bins between only 0 and 120 Hz?

What I tried:

I tried doing something like this by using np.linspace to specify the frequencies (instead of fft.fftfreq), but this introduced some other bug, namely the spectrograms always looked like a mirror, with high power at the highest and lowest frequencies and low power in the middle of the graph, regardless of the range. I'm honestly not sure why this happens with linspace and not with fftfreq. They both return arrays of floats.

Code sample included below. Thanks for any help in advance! Cheers

X = samples_filt  # long 1-D vector of low-pass filtered samples
fs = 1000
time = .05  # window length in sec
N = int(fs * time)  # num samples

tot_len = 5 * fs  # 5 sec of the whole signal
X = X[:tot_len]

f = fftpack.fftfreq(N, 1.0/fs)
# f = np.ceil(np.linspace(0, max_freq, 52)[1:51])  # introduced error described above
mask = (f > 0)  # mask for positive freqs

n_max = int(np.ceil(X.shape[0] / N))  # the number of segments of length N in the sample array data
f_values = np.sum(1 * mask)  # how many values meet mask reqs
spectogram_data = np.zeros((f_values, n_max))

window = sp.signal.blackman(N)  # taper used to improve contrast of spectrogram

for n in range(0, n_max):
    subdata = X[(N * n):(N * (n + 1))]
    F = fftpack.fft(subdata * window)
    spectogram_data[:, n] = np.log(abs(F[mask]))

fig, ax = plt.subplots(1, 1, figsize=(8, 6))
p = ax.imshow(spectogram_data, origin='lower', 
    extent=(0, X.shape[0] / fs, 0, max(f)), 
    aspect='auto',
    cmap=mpl.cm.RdBu_r)
cb = fig.colorbar(p, ax=ax)

Sample output:

For instance, it would be great to have this same graph with 24 frequency bins going up to 120 Hz.

Answer 1

There is a fundamental limit to how much linear frequency resolution you are going to get for a given window size, at least using traditional linear techniques and making no assumptions about the source signal.

There are techniques that does pretty much the same as an FFT but returning only a subset of frequency bins. They are typically used for cost reasons.

For convenience, you might consider resampling your signal to 2x120=240Hz, then do a spectrogram of longer than 50ms windows depending on the temporal resolution you need.

Frequency resolution is: d_f = fs/L

For d_s = 120/50=2.4 Hz you might want something like:

L = fs/d_f=240/2.4=100 samples

Ie a window length of: 100/240 = 0.42 seconds

Answer 2

How could I, for instance, do the FFT specifying something like 50 frequency bins between only 0 and 120 Hz?

If you want a resolution of $\Delta f = \frac{120\mathrm{Hz}}{50}$, then you need sample lengths of at least $\Delta t = \frac{1}{\Delta f} \simeq 0.41\bar6\mathrm s$. With $\Delta t = 50 \mathrm{ms}$ windows, you're limited to $\Delta f = \frac{1}{50\mathrm{ms}} = 20 \mathrm{Hz}$

Recording at 1kHz and throwing away 3/4 of your results to filtering is actually not that bad. You could mess around with pre-filtering, but you'd probably end up finding out that doing so involves FFTs, IFFTs and overlap-adds -- so you may as well just cut out the middle man and do an FFT on the "whole" data.

If you're concerned about getting all the extra stuff from your FFT, there are flavors of the FFT algorithm that are designed to take a vector of real (not complex) samples and return just the positive-frequency side of the FFT. They run a bit faster, because with real input data there are some symmetries that can be leveraged.

Any decent FFT package (i.e. FFTW) has them -- you just need to dig in the documentation a bit to find them.

Answer 3

This is doable with a rowwise implementation of STFT, as opposed to the standard columnwise. It is currently TODO in ssqueezepy, with an already working variant asserting equal against the standard.

Note, however, simply more frequency bins doesn't necessarily yield greater freq resolution. Same reasoning as this answer, but along frequency axis - plus one imposed fundamentally by the resolution of plain DFT.