%%% discoTX.m --- 
%% 
%% Filename: discoTX.m
%% Description: Matlab/octave model illustrating the discontinuous 
%               symbol-synchronized TX processing
%% Author: Juha Yli-Kaakinen
%% Maintainer: 
%% Created: Wed Aug 14 13:24:52 2019
%% Version: 1.0
%% Last-Updated: Tue Nov 10 17:31:52 2020
%%           By: Juha Yli-Kaakinen
%%     Update #: 1553
%% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function discoTX

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

vec = @(x) x(:); % Vectorization function

ran_seed = 10;
if isOctave
  rand('seed', ran_seed);
  randn('seed', ran_seed);
  fontSize = 14;

else
  rng(ran_seed); 
  fontSize = 12;
 
end
set(0, "defaultaxesfontsize", fontSize)
set(0, "defaulttextfontsize", fontSize)
hextodec = @(x) hex2dec(reshape(x,2,3).').'/255;

i0 = 1;                       % Index of the 1st element of the vector

figure(1), clf

Lact = 1*12;                  % Number of active subcarriers (one PRB)
bOFDM = 14;                   % Number of OFDM symbols 

% Input data on active subcarries. Shift around DC
Xact = 1/sqrt(2)*(sign(randn(Lact, bOFDM)) + 1j*sign(randn(Lact, bOFDM)));  

Bin0 = 113;                   % Center bin of the subband

N = 1024;                     % Long FC SFB transform (IFFT) size
  
OLF = 1/2                     % Overlap factor in FC processing
SCS = 15e3;                   % OFDM subcarrier spacing (SCS) in Hz
BS_FC = SCS;                  % FC processing bin spacing in Hz
  

assert(OLF == 1/2, 'In this model, the overlap factor has to be equal to half')

for iModel = [1:3]
  switch iModel
   case 1
    % Continuous FC-OFDM TX
    L = 128; % Short FC transform (FFT) size

    lineColor = hextodec('0072bd'); lineStyle = 'o-';
    legStr = ['Continuous FC-F-OFDM TX, L = L\_OFDM = ' int2str(L)];

   case 2 
    % Discontinuous FC-OFDM TX
    L = 128; % Short FC transform (FFT) size
    
    lineColor = hextodec('77ac30'); lineStyle = 'o-';
    legStr = ['Discontinuous FC-F-OFDM TX, L = L\_OFDM = ' int2str(L)];
    
   case 3 
    % Extrapolating FC-OFDM TX
    
    % Short FC transform (FFT) size. Here 12-bin transition bands are used on the 
    % both sides of the active subcarries, therefore, the minimum L is 36.
    L = 24; 
    
    lineColor = hextodec('a2142f'); lineStyle = 'x-';
    legStr = ['Extrapolating FC-F-OFDM TX, L = L\_OFDM = ' int2str(L)];

   otherwise
    error('Unsupported')
    
  end
  
  I = N/L;                    % Interpolation factor [Eq. (4)]
  L_OFDM = L;                 % Low-rate OFDM transform size
  N_OFDM = I*L_OFDM;          % High-rate OFDM transform size
  Fs = SCS*N_OFDM;            % Output sample rate in Hz
  
  % Determine the FC processing block partition sizes [based on Eqs. (2) and (3)]
  [Ls, Ll, Ns, Nl] = FCpartitioning(L, N, OLF);
 
  % Determine the CP lengths (resempling the LTE and 5G-NR CP lengths for 15 kHz SCS)
  [L_CP, N_CP] = CPlenghts(L_OFDM, N_OFDM, bOFDM);

  disp([legStr ':'])
  disp(['High-rate IFFT, N_OFDM = [' int2str(N_OFDM.') ']'])
  disp([' Low-rate IFFT, L_OFDM = [' int2str(L_OFDM.') ']'])
  disp(['    High-rate CP, N_CP = [' sprintf('%2i  ', N_CP) ']'])
  disp(['     Low-rate CP, L_CP = [' sprintf('%2i  ', L_CP) ']'])
  
  % Number of samples to be processed 
  nSamp_lowrate  = sum(L_CP, 1) + L_OFDM.*bOFDM; % Samples on the input rate
  nSamp_highrate = sum(N_CP, 1) + N_OFDM.*bOFDM; % Samples on the output rate
  
  % Map Lact active subcarries to L_OFDM bins of the OFDM IFFT
  X = circshift([fftshift(Xact, 1); zeros(L_OFDM(1)-Lact, bOFDM)], ...
		[-floor(Lact/2) 0]); 

  % Processing matrices 
  W_N  = 1/sqrt(L)*fft(eye(N)); % DFT matrix of size N x N
  W_L  = 1/sqrt(N)*fft(eye(L)); % DFT matrix of size L x L
  WI_N = sqrt(N)*ifft(eye(N));  % IFFT matrix of size N x N
  WI_L = sqrt(L)*ifft(eye(L));  % IFFT matrix of size L x L
  P    = fftshift(eye(L), 1);   % FFT shift matrix 
  
  % Time-domain analysis window used for overlap-and-add (OSA) processing [Eq. (10a)]
  A_mr = sparse(Ll+(1:Ls), Ll+(1:Ls), ones(1, Ls), L, L);

  % Mapping matrix maps L input bins to N output bins [Eq. (9a)]
  M_LN = @(theta) sparse(mod(ceil(-L/2:L/2-1)+Bin0, N)+1, 1:L, ...
			 theta*ones(1, L), N, L);
  
  % Frequency-domain window function: Ones on the passband and zeros on the stopband 
  % as well as 'nTB' transition-band weights on both sides of the passband
  nTB = 4;                      % Number of transition band weights in FC processing
  excessBW = 2;                 % Passband extension in frequency-domain bins
  win = genWin(nTB, excessBW, Lact, 1, L);
  
  % OFDM modulation, CP inclusion, and zero padding to length 3/2L
  [Z, x_cp] = OFDMmod(X, Lact, bOFDM, L_OFDM, L_CP, L, Ls, Ll);  
  
  switch iModel
   case 1
    %% CONTINUOUS FC PROCESSING (overlap-and-add (OLA) block partitioning)
    
    % Zero-padding needed in FC processing
    prePad = L-Ls; postPad = L-Ls; 
  
    % Total number of samples including the zero padding
    nTotSamp = nSamp_lowrate + prePad + postPad;
    
    % Resulting number of FC prosessing blocks
    nBlks = ceil((nTotSamp-L)./Ls)+1;
    
    % Construction of block diagonal matrix with overlapping blocks
    Fm = zeros(N+Ns*(nBlks-1), L+Ls*(nBlks-1));    

    % Phase rotation [Eq. (9b)]
    theta = @(r) exp(j*2*pi*r*Bin0*Ls/L); 
    
    for r = 0:nBlks-1 % Do for all processing blocks
      % Block-wise OLA processing matrix of size Ns x L [Eq. (7a)]
      Fmr = WI_N*M_LN(theta(r))*diag(win)*P*W_L*A_mr;
      
      % Block-diagonal matrix with overlaping blocks [Eq. (5b)]
      Fm(r*Ns+(1:N), r*Ls+(1:L)) = Fm(r*Ns+(1:N), r*Ls+(1:L)) + Fmr;
      
    end % END for r
    
    % Zero-padding and FC processing [Eq. (5a)]
    w_m = sqrt(I)*Fm*[zeros(prePad, 1); x_cp; zeros(postPad, 1)];
    
    delay = N/L*prePad;
    
   case {2, 3}   
    %% DISCONTINUOUS SYMBOL-SYNCHRONIZED FC TX PROCESSING

    % Segmentation to blocks of length L [Eq. (12)]
    R_Lsm0 = [eye(L); zeros(L/2, L)];
    R_Lsm1 = [zeros(L/2, L); eye(L)];
    z_0 = R_Lsm0.'*Z; 
    z_1 = R_Lsm1.'*Z;
    
    w_mn = zeros(3/2*N, bOFDM);
    
    % Phase rotation [Eqs. (9b), (14b), and (14c)]
    theta = @(r) exp(j*2*pi*r*Bin0*Ls/L); 
    phi = @(n) exp(j*2*pi*Bin0*sum(N_CP((1:n)+i0))/N);
    
    % Selection matrices 
    for n = 0:bOFDM-1
      
      % Time-domain synthesis matrix for 1st block in OLA processing [Eq. (13)]
      A_m0 = sparse(Ll-L_CP(n+i0)+(1:Ls+L_CP(n+i0)), Ll-L_CP(n+i0)+(1:Ls+L_CP(n+i0)), ones(1, Ls+L_CP(n+i0)), L, L); % Selection matrix
      
      % Block-wise OLA processing matrices of size Ns x L [Eq. (7a)]
      F_m0 = sqrt(I)*WI_N*M_LN(theta(0)*phi(n))*diag(win)*P*W_L*A_m0; 
      F_m1 = sqrt(I)*WI_N*M_LN(theta(1)*phi(n))*diag(win)*P*W_L*A_mr; 
      
      % Filtering [Eq. 14a]
      w_m0 = F_m0*z_0(:, n+i0);
      w_m1 = F_m1*z_1(:, n+i0);
      
      % Overlap-and-add of the filtered halfs [Eq. 15]
      R_Nm0 = [eye(N); zeros(N/2, N)]; 
      R_Nm1 = [zeros(N/2, N); eye(N)];
      w_mn(:, n+i0) = R_Nm0*w_m0 + R_Nm1*w_m1; 
      
    end % END FOR n
    
    % Length of the output vector (Eq. 16d)
    Q = N/2 + bOFDM*N_OFDM + sum(N_CP(2:bOFDM));
    w_m = zeros(Q, 1);

    % Concatenation of the filtered symbols [Eq. (16)]
    for n = 0:bOFDM-1
      % Starting index of the nth symbol [Eq. (16b)]
      sigma_n = n*N_OFDM(1) + sum(N_CP((1:n)+i0));

      % Alignment (Eq. 16c)
      Gamma = [ 
	  zeros(sigma_n, 3/2*N)
	  eye(3/2*N)
	  zeros(Q-3/2*N-sigma_n, 3/2*N)
	      ];
      
      % Overlap-and-save (Eq. 16a) 
      w_m = w_m + Gamma*w_mn(:, n+i0);
    
    end % END FOR n
    
    delay = Nl-N_CP(1);
    
  end % END SWITCH

  % Active subcarries on TX side
  iAct_TX = mod([1:Lact]-floor(Lact/2)-1, L_OFDM)+1;
    
  %% RX BASIC OFDM PROCESSING  
  % Throw some samples away (due to zero-padding and overlapping)
  w_hat_m = w_m(delay+(1:nSamp_highrate)); 
  
  % Demodulation back to zero frequency 
  w_hat_m = w_hat_m.*exp(-j*2*pi*Bin0.*BS_FC/Fs*((0:nSamp_highrate-1)+delay).');
  
  % OFDM demodulation                                   
  Y = OFDMdem(w_hat_m, bOFDM, N_OFDM, N_CP);
  
  % Active subcarries on RX sides
  iAct_RX = mod([1:Lact]-floor(Lact/2)-1, N_OFDM)+1;  

  % Difference between the transmitted and received symbols
  E = Y(iAct_RX, :) - X(iAct_TX, :);

  % Power spectral density estimate
  subplot(121)  
  nGridPoints = 2^15;
  if iModel == 5
    [Y, freq] = PSDeval(x_cp, Fs/I, nSamp_lowrate, nGridPoints);
    h = plot(freq/1e6, 10*log10(Y), 'Color', 'm', 'DisplayName', 'Plain CP-OFDM');
    hold on
   
  end
  [Y, freq] = PSDeval(w_m, Fs, nSamp_highrate, nGridPoints);
  
  h = plot(freq/1e6, 10*log10(Y), 'Color', lineColor, 'DisplayName', legStr);
  hold on
  
  % Mean-squared error on active subcarries
  subplot(122)
  h = plot(10*log10(MSE(E)), lineStyle, 'Color', lineColor, 'DisplayName', legStr);
  set(h, 'MarkerFaceColor', 'w', 'lineWidth', 1.4, 'MarkerSize', 8)
  hold on
  
end % END for iModel

subplot(121)
title({['PSD of the TX signal, center frequency = ' int2str(SCS*Bin0/1e3) ' kHz']})
xlabel('Frequency in MHz'), ylabel('PSD in dB')
hold off
legend('Location', 'South')
ylim([-100 10])
xlim([-Fs Fs]/2/1e6)
grid on

subplot(122)
legend('Location', 'North')
hold off
ylabel('MSE error in dB'); xlabel('Active subcarrier n'), grid on
title('MSE error on active subcarriers')
grid on

end

% SUBFUNCTIONS
% ------------------------------------------------------------------------------------
function mse = MSE(Err)

  % Mean-squared error over active subcarriers

  mse = sum(abs(Err).^2, 2)/size(Err, 2);

end

% ------------------------------------------------------------------------------------
function [X_CP_OFDM, x_CP_OFDM] = OFDMmod(X, Lact, bOFDM, L_OFDM, L_CP, L, Ls, Ll)

  % OFDM modulation, CP inclusion, and zero padding to length 3/2xL

  % CP-OFDM symbol vector (used for basic FC processing and PSD evaluation)
  x_CP_OFDM = zeros(bOFDM*L_OFDM+sum(L_CP), 1); 

  % CP-OFDM symbol matrix of size (3/2xL x bOFDM) 
  X_CP_OFDM = zeros(3/2*L, bOFDM);
  
  insCP = @(N, K) ([zeros(K, N-K) eye(K); eye(N)]); % CP insertion matrix
    
  for iSymb = 1:bOFDM
    % OFDM modulation [Eq. (1a)]
    x_OFDMn = sqrt(L_OFDM)*ifft(X(:, iSymb));

    % Add cyclic prefix [Eq. (1b)]
    x_CP_OFDMn = insCP(L_OFDM, L_CP(iSymb))*x_OFDMn;
    
    % Indices of the CP-OFDM symbol
    inx = (iSymb-1)*L_OFDM + sum(L_CP(1:iSymb-1)) + (1:L_OFDM+L_CP(iSymb)); 

    x_CP_OFDM(inx) = x_CP_OFDMn;
    
    % Zero-padding to length 3/2xL (Eq. 10)
    S_L = Ll-L_CP(iSymb);
    S_T = L-Ll-Ls;
    ZP = [ 
	zeros(S_L, L_OFDM+L_CP(iSymb));
	eye(L_OFDM+L_CP(iSymb));
	zeros(S_T, L_OFDM+L_CP(iSymb));
	 ];
    X_CP_OFDM(:, iSymb) = ZP*x_CP_OFDMn;
      
  end % END for iSymb

end % END for iModel

% ------------------------------------------------------------------------------------
function Y = OFDMdem(z, bOFDM, L_OFDM, L_CP)

  % Removal of the CP and OFDM demodulation

  Y = zeros(L_OFDM, bOFDM); 
  for iSymb = 1:bOFDM
    % Here, the reference symbol timing is in the middle of the CP
    time_offset = L_CP(iSymb)/2;
    
    OFDMsymb = z(L_OFDM*(iSymb-1)+sum(L_CP(1:iSymb))-time_offset+(1:L_OFDM));
    
    Y(:, iSymb) = 1/sqrt(L_OFDM)*fft(circshift(OFDMsymb, [-time_offset 0]));    
  
  end % END for iSymb

end

% ------------------------------------------------------------------------------------
function win = genWin(nTB, excessBW, Lact, SCS_fac, L)

  % Frequency-domain window function
  
  % Optimized transition-band weights. Here, the optimization 
  % is carried out following the techniques described in [2]

  tap = [ % JYK; 27.11.19, 14:25
              -0.011485978070882
              +0.238202517933228
              +0.827879691252344
              +1.010936145005123
	];

  win = [ones(ceil(Lact/2*SCS_fac+excessBW), 1);  
	 tap(nTB:-1:1); 
	 zeros(L-Lact*SCS_fac-2*excessBW-2*nTB, 1); 
	 tap(1:nTB); 
	 ones(floor(Lact/2*SCS_fac+excessBW), 1)];

  win = fftshift(win);
  
end

% ------------------------------------------------------------------------------------
function [Y, freq] = PSDeval(y, Fs, nSamp, nGridPoints)

  % Evaluate power spectral density estimate (magnitude-squared response)

  Y = 1/sum(nSamp)*abs(fftshift(fft(y, nGridPoints))).^2;

  freq = ((0:nGridPoints-1)/nGridPoints-0.5)*Fs;
  
end

% ------------------------------------------------------------------------------------
function [L_CP, N_CP] = CPlenghts(L_OFDM, N_OFDM, bOFDM)

  % CP lengths on low- and high-rate sides

  L_CP = repmat(floor(9/128*L_OFDM), [bOFDM 1]); % Low-rate CP lengths
  N_CP = repmat(floor(9/128*N_OFDM), [bOFDM 1]); % High-rate CP lengths
  L_CP(1:7:end) = floor(10/128*L_OFDM);          % Longer CP for 1st symbols in a slot
  N_CP(1:7:end) = floor(10/128*N_OFDM);          % Longer CP for 1st symbols in a slot

end
% ------------------------------------------------------------------------------------
function [Ls, Ll, Ns, Nl] = FCpartitioning(L, N, OLF)

  % FC processing block partition sizes

  Ls = round(L*(1 - OLF)); % Non-overlapping samples for short FC processing transf.
  Ll = ceil((L - Ls)/2);   % Leading overlapping samples for short FC process. transf.
  Ns = round(N*(1 - OLF)); % Non-overlapping samples for long FC processing transf.
  Nl = ceil((N - Ns)/2);   % Leading overlapping samples for long FC process. transf.

end
% ------------------------------------------------------------------------------------
