Combine FFT's of shorter length than sample data to get spectrum of all data

Question

I have a set of sample data, a logarithmic swept sine as it happens, that is 1 second in duration sampled at 100kHz giving 100k samples. I want to calculate the frequency spectrum of the complete data however the FFT library function available only supports 4096 samples maximum.

Is it valid to sequentially FFT 4096 samples and sum the resulting data. Thereby ending up with 4096/2 frequency points that show the frequency content of the entire 100k samples?

I guess it is necessary to window each block of 4096 samples to be FFT'd but this would mean that the portions of the frequency content would be reduce in amplitude artificially.

How does one handle the phase information?

Answer 1

Yes you can do it, but you will loose spectral resolution. The method you mention is a variant of the Welch's spectogram estimation algorithm which relies on the block based estimation of power spectrum of a long sequence in terms of shorter blocks within it.

Given that you have a data of $100$k in length which was sampled at $F_s = 100$k, then a full length DFT after rectangular windowing would provide you a spectral resolution of about $1$ Hz to $3$ Hz, where as a window of $4096$ samples would yield about $25$ Hz to $75$ Hz of spectral resolution.

If this resolution is enough for your purposes, then you can (and indeed you should better) perform the block based averaged spectral analysis in order to trade-off between the reduced spectral resolution and an increased estimation accuracy (decreased estimation variance)

Answer 2

There is a section in Rabiener and Gold’s 1975 DSP book listed near the index term “twiddle factors” around page 347 that shows how to decompose a DFT in 2 dimensions. I used to use it a lot prior to FFTW to do oddball DFT sizes like 384 without zero padding.

Theory and application of digital signal processing
Prentice-Hall signal processing series Lawrence R. Rabiner, Bernard Gold, Prentice-Hall, 1975

First you need to factor the desired DFT size into 2 numbers M and N and then allocate a M by N matrix. Write you data to the matrix in row order skipping to the next row as you fill. This is an in place algorithm.

Do M point DFTs on each column. Next Multiply each element in the matrix by its twiddle factor $exp( j 2 \pi (m-1)(n-1))$ where m and m are the index terms given the Fortran convention. Do N point DFTs on each row. Read out the results in column order.

Matlab code:

clear all
M=3;
N=32;
x=linspace(0,10,M*N);
X=reshape(x,N,M).'; % read in as rows
Twiddle=zeros(size(X));
for i=1:M
for k=1:N
Twiddle(i,k)=exp(-1j*2*pi*(i-1)(k-1)/(NM));
end
end
X=fft(X) % fft on each column
X=X.*Twiddle;% element by element product
X=fft(X.').' ; %fft on each row
y=reshape(X,N*M,1); % read out as columns
figure(1)
plot(abs(y),'linewidth',2)
title('Composite DFT')
figure(2)
plot(abs(fft(x)),'linewidth',2)
title('Direct DFT')
figure(3)
plot(x,'linewidth',2)
title('time series')

Plots:

There is also another order mentioned in Rabiener’s book. Twiddle first, row fft, column fft. He also makes it a point to note that each row or column dft can be implemented, in turn as a composite square as long as N is not a prime number. The row in, column out indexing is also interesting because it generalizes bit reverse addressing.

In summary, you can do a 100k DFT with 4k DFTs if you have enough memory.