Computation of Only Even or Odd Frequency Bins of DFT

Question

I have an algorithm where I am computing the FFT of a large signal. However, I desire only the even or odd terms of the DFT of the signal, but not both. Currently, I discard these undesired terms. Is there a way by which I may appreciably reduce the runtime complexity when only every other DFT term is desired?

How might I begin to find such an optimization? If it is unclear, I will look for how the larger, encompassing algorithm could be improved, but I'm starting by examining the smallest, most expensive portion before thinking about more fundamental changes to the system.

Answer 1

Consider a sequence $x[n]$ of length N, and assume $X[k]$ is its N-point DFT given by $$X[k] = \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N}nk}}$$ $k = 0, 1,..., N-1$, which is computed through an N-point FFT. If you wish to compute the even-indexed ($k=0,2,4,...$), or the odd-indexed ($k=1,3,5,...$) indexed samples of $X[k]$, you can proceed with the following:

Denote the even indexed samples of $X[k]$ as $X_e[k]$ : \begin{align} X_e[k] = X[2k] &= \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N}n(2k)}} &\scriptstyle{\text{term 2k computes even samples of X[k]}}\\ X_e[k] &= \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{factor 2 is moved into N/2 term}}\\ X_e[k] &= \sum_{n=0}^{\frac{N}{2}-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} + \sum_{n=\frac{N}{2}}^{N-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{Divide the sum into 2 halves}}\\ X_e[k] &= \sum_{n=0}^{\frac{N}{2}-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} + \sum_{n=0}^{\frac{N}{2}-1}{x[n+\frac{N}{2}]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{adjust the 2nd sum range}}\\ X_e[k] &= \sum_{n=0}^{\frac{N}{2}-1}{ (x[n]+x[n+\frac{N}{2}])e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{merge the 2 sums}}\\ \end{align}

Final form suggest that the required even samples of N-point DFT of the sequence $x[n]$ can be computed from an N/2-point DFT of a new sequence $x_{half}[n]= x[n]+x[n+ N/2]$ whose length is half that of $x[n]$, assuming N is even.

The signal $x_{half}[n]$ is simply computed by adding the first half of the signal $x[n]$ into its second half. (if the length N of signal x[n] is not even padd a zero to its tail to make it even)

The following matlab/octave code gives you the desired even indexed samples of X[k] without computing a full N-point DFT of x[n]:

Xe = fft ( x(1:N/2) + x(N/2 + 1 : end), N/2);

Note that because of the addition of the halves before the FFT, efficiency will degrade from a pure N/2-point FFT.

The case for the odd indexed samples proceeds similary but results in a more complex form which I would like to summarize in the following Matlab script:

% Let x[n] be the signal of length N
xc = x .* exp(-j*2*pi*[0:N-1]/N) ; % Multiply x[n] with a complex phase term
Xo = fft (xc(1:N/2) + xc(N/2 + 1, N), N/2); % odd indexed samples of X[k]

Due to a complex multiplication before the FFT, in this case of computing odd indices performance is further reduced but still preferable over a direct N-point FFT.

The following Matlab/Octave excerpt demonstrates the computation of both the even and the odd indexed samples.

N = 1024;           % original sequence length
x = randn(1,N);
X = fft(x,N);       % original DFT of X

Xe = fft( x(1:N/2) + x(N/2 + 1: N) , N/2);                          % Even   samples only from N/2 point FFT
Xo = fft( ( x(1:N/2) - x(N/2+1:N) ).*exp(-j*pi*[0:2:N-2]/N) , N/2); % Odd samples only from N/2 point FFT

figure,stem( abs ( Xe - X(1:2:N)) ); title('Xe - X(1:2:N)');
figure,stem( abs ( Xo - X(2:2:N)) ); title('Xo - X(2:2:N)');

Due to numerical roundoff effects, the difference plotted by the stems will be nonzero but very small...