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 !
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;
- FFT of Pure Real Sequences, Engineering Productivity Tools Ltd.
4 thoughts on “Computing DFT of real signals in MATLAB”
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)
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.