%%% FFTdecomp.m --- 
%% 
%% Filename: FFTdecomp.m
%% Description: Simplified discontinuous processing for L = 16 as described by (25) in [1]. 
%% Author: Juha Yli-Kaakinen
%% Created: Tue Jan  7 14:37:48 2020
%% Version: 1.0
%% Last-Updated: Fri Jun 12 13:15:24 2020
%%           By: ylikaaki
%%     Update #: 461
%% Keywords:
%% References:
%% [1] Juha Yli-Kaakinen, AlaaEddin Loulou, Toni Levanen, Kari Pajukoski, Arto Palin, 
%%     Markku Renfors, and Mikko Valkama, "Frequency-domain signal processing for 
%%     spectrally-enhanced CP-OFDM waveforms in 5G new radio," submitted IEEE
%%     Transactions on Wireless Communications, 2020.
%% [2] H. Sorensen, M. Heideman, and C. Burrus, "On computing the split-radix FFT,"
%%     IEEE Trans. Acoust., Speech, Signal Processing, vol. 34, no. 1, pp. 152-156, 
%%     Feb. 1986.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all

L = 16;

W = @(L) fft(eye(L)); WI = @(L) ifft(eye(L));   % DFT and IDFT matrices

isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; 

ran_seed = 10;
if isOctave
  rand('seed', ran_seed);
  randn('seed', ran_seed);
  
else
  rng(ran_seed); 
 
end
% Number of real multiplications for split-radix FFT of length N [2]
multSR = @(N) N*(log2(N) - 3) + 4;

% Generate some random data of 3/2*L samples corresponding nth OFDM symbol 
% (i.e., two blocks of L samples with 50 % overlap)
gm = randn(3/2*L,1) + j*randn(3/2*L,1);

% Overlapping blocks of length L in two vectors
gm0 = gm(1:L);  gm1 = gm(L/2+(1:L));

%% Original implementation having two IFFTs and one FFT of length L
S_Ls = [zeros(L/2, L/4) eye(L/2) zeros(L/2, L/4)];
x_mn_0 = fft([S_Ls*ifft(gm0); S_Ls*ifft(gm1)]);
cmplx_0 = 3*multSR(L); % Corresponding complexity
disp(['No. real multiplications for the original approach is ' int2str(cmplx_0)])

%% Implementation according to Eq. (23)
S = [zeros(L/4, L); zeros(L/2, L/4) eye(L/2) zeros(L/2, L/4); zeros(L/4, L)]; % Left-hand of Eq. (19) 
P = @(L, shift) circshift(eye(L), [0 shift]); % Eq. (8)
x_mn_1 = W(L)*[P(L, L/4)*S*WI(L)*gm0 + P(L, -L/4)*S*WI(L)*gm1]; % Eq. (23)
cmplx_1 = 3*multSR(L);
disp(['No. real multiplications for implementation in Eq. (23) ' int2str(cmplx_1)])

assert(all(abs(x_mn_0-x_mn_1) < 1e-10))

%% Implementation according to Eq. (24)
F_0 = W(L)*circshift(eye(L), [0 +L/4])*S*WI(L); % Left-hand subequation of Eq. (24b)
F_1 = W(L)*circshift(eye(L), [0 -L/4])*S*WI(L); % Right-hand subequation of Eq. (24b)
x_mn_2 = F_0*gm0 + F_1*gm1; % Eq. (24a)
cmplx_2 = 4*multSR(L);
disp(['No. real multiplications for implementation in Eq. (24) ' int2str(cmplx_2)])

assert(all(abs(x_mn_0-x_mn_2) < 1e-10))

%% Implementation according to Eq. (25) (one FFT and one IFFT of length L)
w = exp(-2i*pi/L);
Omega = @(r) diag(w.^vec(r*(0:L-1))); % Eq. (25c)
S_hat = P(L, L/4)*S*P(L, -L/4); % Eq. (25d)
F_sub = W(L)*S_hat*WI(L); F = Omega(L/4)*F_sub*Omega(-L/4); % Eq. (25b)
x_mn_3 = Omega(-L/4)*[F*(gm0-Omega(L/2)*gm1) + Omega(L/2)*gm1]; % Eq. (25a)
cmplx_3 = 2*multSR(L);
disp(['No. real multiplications for implementation in Eq. (25) ' int2str(cmplx_3)])

assert(all(abs(x_mn_0-x_mn_3) < 1e-10))



