If compute time doesn't matter, you can overcome the memory problem easily with a dual for-loop implementation. But there's an FFT alternative with caveats: break up signal into M seqments, take M FFT's, add them up. The result is exact, limited in the lowest possible frequency it can represent. You could, however, do both; FFT for high frequencies, and for-loop for leftovers.
FFT in segments
Recall that FFT computes DFT, which in turn multiplies input by sines of cosines of varying frequencies. Let's take x of length 512, f=4 for left half, f=64 for right:
real(FFT(x)[1]), for example, was computed by multiplying:
And imag(FFT(x)[1]) would be computed by replacing the orange cosine of f=1 with orange sine of f=1. Now, what if N=512 is too long, and we decide to break it up into M=2 frames? Then, from those FFT's (FFT256) perspective, real(FFT(x)[1]) looks like
and from the FFT512's perspective, this looks like
which is real(FFT(x)[2]). In other words, FFT(x[:256])[1] + FFT(x[256:])[1] == FFT(x)[2].
Now imagine the highest frequency basis (sine or cosine) in either FFT256 or FFT512; what's the limit? It's [-1, 1, -1, 1], sample-to-sample (k=1/N). No matter how long the FFT is, this basis is always the greatest (normalized) frequency possible from x's point of view. As from above example, shortening segments corresponded to being unable to compute FFT(x)[1] - because that'd require a k=0.5 basis in FFT256:
It's important to understand "normalized frequency" to correctly combine the FFT's later. But basically, numerically the shorter segment computations are exactly the same as longer segment computations with higher-frequency bases.
For a billion point signal, the highest FFT basis frequency is k=500,000,000. What this actually means, is that the dynamic range of frequencies you can capture is 500,000,000:1 (look up normalized frequency). So, if in data you expect the ratio between the highest and lowest frequency to be 10 times less, then you can capture the entire spectrum by splitting up the billion points into 10 segments and adding up the results. If it's any more, you can capture a part of it with the above method, and another part with for-loops.
Additional notes:
- If the input is real-valued, be sure to use
rfft (real FFT), which takes half the memory.
- Beware of any solutions proposing a window of any kind other than rectangular; they will distort the spectrum (but it may be practically insignificant).
- Beware zero-padding
- For very long DFT, be sure to use float64, and even that may not be enough depending on level of accuracy you seek
For-loops rDFT implementation: (Python but easy to port to MATLAB)
def dft(x): # unnormalized
N = len(x)
reX = np.zeros(N // 2 + 1)
imX = np.zeros(N // 2 + 1)
for k in range(N // 2 + 1):
for n in range(N):
reX[k] += x[n] * np.cos(2 * np.pi * k * n / N)
imX[k] += x[n] * np.sin(2 * np.pi * k * n / N)
return reX - imX * 1j
Additionally
If working in Python, I recommend Numba for CPU-parallelism, and Cupy to run code on GPU. Both help with overcoming Python's slug for-loops and Numpy's single-threadedness.
