# FFT for long waveform

I should take FFT of a very long vector of length 1,636,399,763 (i.e. length order is 10^9) preferably in MATLAB. However, on 128 GB RAM, my code throws Out of Memory error on line part2Vec = repmat(part2Vec_OneRepetition, 1000, 1) from script below which is attempting to build a vector of length 1,636,763,000 before even attempting to build totalVec which is of length 1,636,399,763 and take FFT which would be of length nFftLen = 2,147,483,648.

Any thoughts how do do this? For instance, can I take FFT in chunks and if yes, how I can do that?

%part1Vec = this is a "complex" column vector of length "1,636,763". Also, "complex" means each element is a complex number a+b*j.
%part2Vec_OneRepetition =  this is a complex column vector of length "1,636,763"
part2Vec = repmat(part2Vec_OneRepetition, 1000, 1); % This is a complex column vector of length "1,636,763,000".
totalVec = [part1Vec ; part2Vec]; % This is a complex column vector of length "1,638,399,763".
nLen = length(totalVec);% This is "1,638,399,763".
nFftLen = 2.^(nextpow2(nLen));
myFft = fftshift(fft(totalVec, nFftLen));

• Rather than figuring out how to perform a 1.6 billion point FFT you should ask yourself if you really need to perform such an FFT. In any case I think you will have trouble. You're expecting built-in trig function software approximations to compute accurate sines and cosines of angles as small as 2 times ten to the minus seven degrees. Can any trig function software do that accurately? Oct 2, 2020 at 10:13
• You must have information beyond the signal length that you can share with us. Oct 6, 2020 at 21:33

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.

• 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).
• 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