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.

4 thoughts on “Computing DFT of real signals in MATLAB”

  1. There are bugs in your code:

    In fftr exp( ) must read:
    exp(((N/4 – 1):-1:-(N/4 – 1))/(N/2) *pi*1i)

    In ifftr exp( ) must read:
    exp(((N/4):(3*N/4 – 1))/(N/2)*pi*1i)

    Regards,
    Clemens

  2. Dear Clemens,
    Thank you for noting this. Obviously, I have tested the above functions with row input vectors only. You must have used a column vector as the input instead.
    I added a pre-check of input vector orientation, to make functions more robust.

  3. Thank you very much simon! This is very useful! For your information im using these to test out a spectral re-meshing procedure in MATLAB for a Pseudo-spectral DNS code!

  4. I’ve made some minor optimization to the code above:
    Since the size of the FFT blocks is usually constant within the same problem, the result of the exp(k * (pi * 1i / (N/2))) is now cached for re-use.

Leave a Reply

Your email address will not be published. Required fields are marked *