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

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

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

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.

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!

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.