# 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.
%

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).
%

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. Clemens Schlegel says:

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. Gabs says:

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.

This site uses Akismet to reduce spam. Learn how your comment data is processed.