FFT of a $N$-length real sequence via FFT of a N/2-length complex sequence

Question

I would like to optimise FFT algorithm by only using a FFT of size $N/2$ to transform a real sequence of length $N$. There are already many questions and many excellent answers on the matter (here for example).

I've been trying to follow the linked answer to implement the splitting algorithm in Matlab and I checked against the standard FFT Matlab routine. I found a shift of 1 sample in the spectrum (sine input) and a weird peak at a (seems to me) random frequency. I would really appreciate a review of my code or simply a hint to debug the problem.

As a side note, I'm trying this route because I need a memory efficient FFT for an embedded device.

close all;
clear all;
clc;

%% Parameters
sample_rate = 50; % Hz
frequency = 10; % Hz
N = 256;
p = 1:N;
t = (p-1) / sample_rate;

f = sample_rate*(0:(N/2 -1))/N; % frequency vector


%% Generate real sequence
data = sin(2*pi*p*frequency/sample_rate);

figure();
plot(t, data);
title('Input');

%% Split even-odd
% Even indeces are odd because 
% Matlab starts indexing from 1
idx_even = 1:2:N;   
idx_odd = idx_even+1;
data_even = data(idx_even);
data_odd = data(idx_odd);
data_cpx = complex(data_even, data_odd);

%% Perform FFT of length N/2
X = fft(data_cpx, N/2);


%% Perform Butterfly
k = 0:(N/2 -1);
W = exp(-2*pi*1i*k/N);

X_even =     0.5*(X(k+1) + conj(X(N/2 - k)));
X_odd  = -1i*0.5*(X(k+1) - conj(X(N/2 - k)));
G = X_even + W .* X_odd;

%% Perform FFT of length N to check results
X_N = fft(data, N);

%% Plot
figure();
plot(f, real(X_N(1:128)));
hold on
plot(f, real(G));
title('Real');

figure();
plot(f, imag(X_N(1:128)));
hold on
plot(f, imag(G));
title('Imag');

figure();
plot(f, abs(X_N(1:128)));
hold on
plot(f, abs(G));
title('Mag');

figure();
plot(f, angle(X_N(1:128)));
hold on
plot(f, angle(G));
title('Phase');

Answer 1

The following Matlab Excerpt corrects your mistake. You had a problem with the DFT length N of the composed sequence... Note that I used a modulus operator for simplicity, you should replace it with indexing for efficiency.

%% Perform Butterfly

%k = 0:(N/2 -1);
%W = exp(-2*pi*1i*k/N);
%X_even =     0.5*(X(k+1) + conj(X(N/2 - k )));
%X_odd  = -1i*0.5*(X(k+1) - conj(X(N/2 - k )));
%G = X_even + W .* X_odd;

k = 0:N-1;
W = exp(-2*pi*1i*k/N);
X_even =     0.5*(X(mod(k,N/2)+1) + conj(X( mod(- k,N/2)+1 )));
X_odd  = -1i*0.5*(X(mod(k,N/2)+1) - conj(X( mod(- k,N/2)+1 )));
G = X_even + W .* X_odd;

fN = sample_rate*(0:N-1)/N;   % new frequency vector data length N...

Mathematically speaking, the line where you made the implementation error is the Even-Odd combination stage:

$$ X[k] = X_e^{N/2}[k] + e^{-j\frac{2\pi}{N}k } X_o^{N/2}[k] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

Where $X_e[k]$ and $X_o[k]$ are $N/2$ point DFTs of the even and odd indexed samples of the original length $N$ real sequence $x[n]$. And $X[k]$ is the N-point DFT of $x[n]$. Note that frequency index $k$ ranges through $0,1,..,N-1$. But those $X_e[k]$ and $X_o[k]$ are half length (but periodic!) sequences. Hence for $k> N/2-1$ they repeat. Remeber that DFT seqeuences are inherently periodic: $X_e[k+N/2] = X_e[k]$. In order to use the base length of those half-length periodic sequences, you can utilize the modulus operator:

$$ X[k] = X_e^{N/2}[(k)_{N/2}] + e^{-j\frac{2\pi}{N}k } X_o^{N/2}[(k)_{N/2}] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

where $(k)_{N/2} = \mod(k , N/2)$ is the modulus of $k$ wrt $N/2$. Argument of every DFT sequence should circularly wrap around...