Computing DFT of real signals in MATLAB

The discrete Fourier transform (DFT) calculates a spectrum of a complex input vector. However, the input vector is usually real, representing a sampled real signal. For real vectors (the imaginary part of all vector components is zero) the upper half of the resulting DFT spectrum is a conjugate mirror of the lower half. Could therefore the computation of the DFT of the real vectors be made with almost half the effort? Fortunately, yes [1]!

Since MATLAB does not provide a function for computation of DFT spectra of real signals out of the box, I developed my own to save time and space in my computations.

fftr.m – Fast Fourier Transform of real vectors

function X = fftr(x)
%FFTR Discrete Fourier transform of real vectors.
%   FFTR(X) is the discrete Fourier transform (DFT) of real vector X.
%
%   For length N input vector x, the DFT is a length N/2+1 vector X, with
%   elements:
%                    N
%      X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                   n=1
%
%   The upper half of discrete Fourier transform of real vectors is
%   conjugate mirror of the lower half and is not computed and returned by
%   this function to save time and space.
%
%   See also FFT, IFFTR.

%   Copyright 2012 Simon Rozman.

persistent E;

% Determine size and orientation of the input vector.
[W, H] = size(x);
N = W * H;
if W == 1
    k =  (N/4 - 1):-1:-(N/4 - 1)  ;
elseif H == 1
    k = ((N/4 - 1):-1:-(N/4 - 1))';
else
    error('Input parameter x is not a vector.');
end
if ~isreal(x)
    error('Input vector is not real.');
end
if any(size(E) ~= size(k))
    E = exp(k * (pi * 1i / (N/2)));
end

% Interleave even and odd samples of real input vector to form complex
% vector half the size, and compute it's DFT using FFT.
X = fft(x(1:2:end) + 1i * x(2:2:end));

% Use DFT of the interleaved complex vector to determine the first N/2+1
% elements (lower half) of the spectrum as if it was computed using full
% length DFT of the real input vector.
Xi         = conj(X(end:-1:2));
X(2:end)   = 0.5 * ( ...
    (X(2:end) + Xi) - ...
    (X(2:end) - Xi) .* E );
X(N/2 + 1) = real(X(1)) - imag(X(1));
X(1)       = real(X(1)) + imag(X(1));

ifftr.m – Inverse fast Fourier Transform of real vectors spectra

function x = ifftr(X)
%IFFT Inverse discrete Fourier transform of real vectors.
%   IFFT(X) is the inverse discrete Fourier transform of X, representing an
%   upper N/2+1 spectral samples of a real vector. Thus, the result is a
%   real vector of length N = 2*(length(X)-1).
%
%   See also IFFT, FFTR.

%   Copyright 2012 Simon Rozman.

persistent E;

% Determine size and orientation of the input vector.
[W, H] = size(X);
if W == 1
    N = 2 * (H - 1);
    k =  (N/4):(3*N/4 - 1)  ;
elseif H == 1
    N = 2 * (W - 1);
    k = ((N/4):(3*N/4 - 1))';
else
    error('Input parameter X is not a vector.');
end
if any(size(E) ~= size(k))
    E = exp(k * (pi * 1i / (N/2)));
end

% Prepare the spectrum as the interleaved complex vector half the size
% would have, and calculate the IDFT using IFFT.
Xi = conj(X(end:-1:2));
X  = 0.5 * ( ...
    (X(1:(end - 1)) + Xi) + ...
    (X(1:(end - 1)) - Xi) .* E );
x  = ifft(X);

% Deinterleave: real values of the interleaved complex vector represent
% even samples of output real vector, imaginary values represent odd
% samples.
x_even = real(x);
x_odd  = imag(x);
x(2:2:N) = x_odd;
x(1:2:N) = x_even;

References

1. FFT of Pure Real Sequences, Engineering Productivity Tools Ltd.