function Xs = fftreal(x)
% This function calculates the FFT of an N-point real signal x[n] by use
% of function fft241.m. First, we break down the given sequence into
% two N/2 point sequences, then use fft241.m and finally the radix-2
% Decimation In Time relations to compute the FFT of the whole sequence.
% Caution: Works well only for N even.
N = length(x);
k = 0:N/2-1;
W = exp(-1i*2*pi*k/N);
x1 = zeros(1,N/2);
y1 = zeros(1,N/2);
even_time_ind = 1:2:N-1; % Form the even sequence x, x, x, ..., x[N-2].
odd_time_ind = 2:2:N; % Form the odd sequence x, x, x, ..., x[N-1].
% Form two REAL sequences x1 and y1:
x1 = x(even_time_ind);
y1 = x(odd_time_ind);
% Use fft241.m to compute the two N/2 DFT's:
[X Y] = fft241(x1,y1);
% Use radix-2, Decimation-In-Time (DIT) relations to compute the DFT of x[n].
Xs = zeros(1,N);
Xs(1:N/2) = X + W.*Y;
Xs(N/2+1:N) = X - W.*Y;