1
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ The scheme that you proposed won't give you the same results as a 100k-point FFT (as Fat32 explains below), but it is possible to build up a larger FFT from multiple smaller FFTs. This sort of decomposition is how the FFT works, actually. $\endgroup$ – Jason R Oct 26 '17 at 15:42
2
$\begingroup$

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)

$\endgroup$
  • $\begingroup$ Noted that the spectral resolution is determined by sample rate and number of samples used for the FFT, that's fine. Now Im just trying to figure out if the phase needs any treatment or is correct by definition. The FFT results in a complex number vextor so is it valid to simply keep adding the resulting FFT complex number vectors as I move through the time domain? $\endgroup$ – happenstance Oct 30 '17 at 0:25
  • $\begingroup$ yes, add complex numbers... $\endgroup$ – Fat32 Oct 30 '17 at 1:05
2
$\begingroup$

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: enter image description here

enter image description here

enter image description here

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.

$\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.