3
$\begingroup$

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');
$\endgroup$
  • $\begingroup$ %% Split even-odd % Even indeces are odd because % Matlab starts indexing from 1 and people wonder why i get so bent outa shape with MATLAB's inability to redefine the array origin from 1 to 0 or any other integer. $\endgroup$ – robert bristow-johnson Oct 10 '18 at 0:51
  • $\begingroup$ @robert bristow-johnson That's pretty annoying, especially when using Matlab only to develop quickly and the real work is in C $\endgroup$ – jack Oct 10 '18 at 6:46
1
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.