3
$\begingroup$

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)); 
$\endgroup$
2
  • 1
    $\begingroup$ 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? $\endgroup$ Oct 2 '20 at 10:13
  • $\begingroup$ You must have information beyond the signal length that you can share with us. $\endgroup$
    – Michael_RW
    Oct 6 '20 at 21:33
3
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Nice. This helps with a problem I’m working on right now. $\endgroup$
    – Ed Tate
    Mar 3 at 1:14
0
$\begingroup$

Perform an overlapped, windowed analysis with e.g. a Hann window and 50% overlap.

Choose an appropriate window size, e.g. 4096 or any other power of two. Since you do not seem to be interested in the spectral content over time, you could opt for high frequency resolution, like 2^15 or something in that range. Then, simply hop over your signal with hops of half the window length, do FFT for each window and add up all FFT results. If you absolute energies matter, divide the final result by the number of windows.

Concerning the actual implementation: I'm not really a matlab guy. If I had to do this, it would be in Python using numpy. A simple generator function for windowing ensures, that at no point in time the whole signal has to be held in memory. I'm sure, something similar exists in matlab.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.