Filtered Multicarrier Transmission

Page

Juha Yli-Kaakinen, Markku Renfors, and Eleftherios Kofidis, and “Filtered Multicarrier Transmission,” Chapter 2, Wiley 5G Ref, John Wiley & Sons, 2020.

Abstract

While orthogonal frequency-division multiplexing (OFDM) has been adopted as the waveform of choice in existing and emerging broadband wireless communication systems, investigations of advanced multicarrier transmission schemes have continued with the aim of eliminating or mitigating the essential limitations of the OFDM scheme. This chapter focuses on multicarrier schemes with enhanced spectrum localization, which manage to reduce the spectral sidelobes, problematic in various advanced communication scenarios. These developments include schemes for enhancing the OFDM waveform characteristics through additional signal processing, as well as filter bank multicarrier (FBMC) waveforms utilizing frequency-selective filterbanks instead of plain (inverse) discrete Fourier transform (DFT) processing for waveform generation and demodulation.

Figure: Channelization of two bandwidth parts with 15 kHz and 30 kHz sub-
carrier spacing (SCS) within a 10 MHz channel. Lower subfigure shows the over-
all responses whereas the upper subfigure shows the transition-band details. The
minimum attenuation at channel edge, As, is given for each waveform

Matlab/octave mode

Matlab/octave model for realizing plain CP-OFDM, time-domain and frequency-domain filtered OFDM processing (TD-F-OFDM and FC-F-OFDM, respectively) as well as windowed overlap-and-add processing (WOLA) is shown below. This model transmits the waveform using one of the above-mentioned TX processing schemes and carries out the RX processing using all the alternatives.

This model can also be downloaded from https://yli-kaakinen.fi/filteredOFDM.m

%% Description: Model for illustrating the channelization of two 5 MHz bandwidth parts 
%%              within a 10 MHz channel using the processing alternatives possible for 5G-NR. 
%% 
%% Author: Juha Yli-Kaakinen (juha.yli-kaakinen@tuni.fi), 
%%         Tampere University, Finland
%% Created: Mon Jan 22 13:11:42 2018
%% Version: $1.0$
%% Last-Updated: Sun Dec  1 15:17:57 2019
%%           By: Juha Yli-Kaakinen
%%     Update #: 4321
%% Keywords: 
%% References:
%% [1] Markku Renfors, Juha Yli-Kaakinen, Toni Levanen, Mikko Valkama, Tero Ihalainen,
%%     and Jaakko Vihriala "Efficient fast-convolution implementation of filtered 
%%     CP-OFDM waveform processing for 5G," in Proc. IEEE GC 2015 Workshop on 5G & 
%      Beyond - Enabling Technologies and Applications (in conjunction with IEEE 
%%     GLOBECOM), San Diego, CA, USA, Dec 6-10, 2015. 
%%     https://doi.org/10.1109/GLOCOMW.2015.7414034
%% [2] Juha Yli-Kaakinen, Toni Levanen, Sami Valkonen, Kari Pajukoski, Juho Pirskanen,
%%     Markku Renfors, and Mikko Valkama, "Efficient fast-convolution based waveform
%%     processing for 5G physical layer," IEEE Journal on Selected Areas in
%%     Communications (Special Issue on Deployment Issues and Performance
%%     Challenges for 5G), vol. 35, no. 6, pp. 1309-1326, 2017. 
%%     https://doi.org/10.1109/JSAC.2017.2687358
%%  Usage: call with filteredOFDM(Model) where Model is either 'CP-OFDM, 'FC-F-OFDM',
%%     'TD-F-OFDM', or 'WOLA'.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function filteredOFDM(TXproc)

if (nargin == 0)
  TXproc = 'FC-F-OFDM';
  
end

ran_seed = 10; 

isOctave = exist('OCTAVE_VERSION', 'builtin') == 5;
if isOctave 
  randn('seed', ran_seed);
else
  rng(ran_seed);
end

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

nAct = [26; 13]*12;         % Number of active subcarriers (25 and 11 physical 
                            %   resource blocks (PRBs) for 1st and 2nd bandwidth
                            %   part (BWP), respectively)
L = 512;                    % Fast-convolution (FC) processing short transform (FFT) 
                            %   size (L_m in [2])
I = 2;                      % Interpolation factor 
N = I*L;                    % FC processing long transform (IFFT) size
L_OFDM_LR = [512; 256];     % TX OFDM IFFT sizes for BWPs giving 15 kHz and 30 kHz 
                            %   subcarrier spacings (SCSs)
Bin0 = 162*[-1; 1];         % Center bins of the BWPs (180 kHz guardband)

Fs = N/L*L_OFDM_LR(1)*15e3; % Output sampling rate (30.72 MHz)
                                                    
bOFDM = [14; 28];           % Number of OFDM symbols (14 for 15 kHz SCS and 28 for 
                            %   30 kHz SCS)
L_OFDM_HR = I*L_OFDM_LR;    % OFDM FFT size on the high-rate side

BS_FC = Fs/N;               % FC processing bin spacing (BS)
SCS_OFDM = Fs/I./L_OFDM_LR; % OFDM processing SCS
SCS_fac = SCS_OFDM/BS_FC;   % Ratio between OFDM SCS and FC processing SCS

L_CP_LR = 9*L_OFDM_LR/128;  % CP length on the low-rate side

% Symbol-wise CP lengths on the low-rate side
L_CP_LR(1, 1:bOFDM(1)) = 9*L_OFDM_LR(1)/128;     % CP length for BWP 1
L_CP_LR(1, 1:7:end)  = L_CP_LR(1, 1:7:end) + 8;  % For the 1st symbols in a slot
L_CP_LR(2, 1:bOFDM(2)) = 9*L_OFDM_LR(2)/128;     % CP length for BWP 2
L_CP_LR(2, 1:14:end) = L_CP_LR(2, 1:14:end) + 8; % For the 1st symbols in a slot

if strcmpi(TXproc, 'UFMC') | strcmpi(TXproc, 'FC-UFMC')
  L_CP_LR = zeros(size(L_CP_LR)); 
  
end

L_CP_HR = I*L_CP_LR;        % CP length on the high-rate side

nTB = 12;                   % Number of transition band weights
excessBW = 3;               % The passband is made wider by 'excessBW' bins

assert(all(nAct < L_OFDM_LR), ...
       'Number active subcarriers has to smaller than OFDM transform size')
assert(all(nAct.*SCS_fac+2*nTB < L), ...
       ['Number of non-zero bins of the FD window has to be smaller than ' ...
        'FC processing short transform'])

OLF = 1/2;                  % Overlap factor 
Ls = round(L*(1 - OLF));    % Non-overlapping part on the low-rate side
Ll = ceil((L - Ls)/2);      % Leading overlapping part on the low-rate side
Ns = round(N*(1 - OLF));    % Non-overlapping part on the high-rate side
Nl = ceil((N - Ns)/2);      % Leading overlapping part on the high-rate side
prePad = L-Ls; postPad = L-Ls;
pre0 = N/L*(prePad-Ll);

% Time-domain overlap-save matrix
S_Ns = sparse(1:Ns, Nl+(1:Ns), ones(1, Ns), Ns, N);
S_Ls = sparse(1:Ls, Ll+(1:Ls), ones(1, Ls), Ls, L);

assert(OLF < 1, 'Overlap factor has to smaller than one')

% Number of samples to be processed 
nSamp = unique(sum(L_CP_LR, 2) + L_OFDM_LR.*bOFDM);

assert(numel(nSamp) == 1, 'Number of samples should be the same for both BWPs')

nTotSamp = nSamp + prePad + postPad;

% Number of FC prosessing blocks
nBlk = ceil((nTotSamp-L)./Ls)+1;

% Preallocate vectors for MSE and PSD evaluations
for iBWP = 1:2
  MSE_fofdm{iBWP} = zeros(nAct(iBWP), 1);
  MSE_ofdm{iBWP}  = zeros(nAct(iBWP), 1);
  MSE_fc{iBWP}    = zeros(nAct(iBWP), 1);
  MSE_wola{iBWP}  = zeros(nAct(iBWP), 1);
  win{iBWP} = zeros(L, 1);
end

BWPs = 1:2;
Err = 0;
nBurst = 100; 


nGridPoints = 2^16;
PSD = zeros(nGridPoints, 1);
for iBWP = BWPs 
  % Frequency-domain (FD) windows used for fast-convolution processing.
  % It should be pointed out that the resolution of the (FD) window 
  % is determined by the FC processing bin spacing
  win{iBWP} = genWin(nTB, excessBW, nAct(iBWP), SCS_fac(iBWP), L);
  
  % Input and output bins of the subband
  iActFC{iBWP}   = mod([1:nAct(iBWP)]-floor(nAct(iBWP)/2)-1, L_OFDM_LR(iBWP))+1;
  iActOFDM{iBWP} = mod([1:nAct(iBWP)]-floor(nAct(iBWP)/2)-1, L_OFDM_HR(iBWP))+1;
  iActUFMC{iBWP} = mod([1:nAct(iBWP)]-floor(nAct(iBWP)/2)-1, nAct(iBWP))+1;
  switch TXproc
   case {'FC-F-OFDM'}
    iActTX{iBWP} = iActFC{iBWP};
    
   case {'CP-OFDM', 'WOLA', 'TD-F-OFDM'}
    iActTX{iBWP} = iActOFDM{iBWP};
    
   otherwise
    error(['model ''' TXproc ''' is not supported'])
    
  end
  
  % Frequency-domain mapping matrix (maps L_m points to N points)
  M_LN{iBWP} = sparse(mod(ceil(-L/2:L/2-1)+Bin0(iBWP), N)+1, 1:L, ones(1, L), N, L);
      
end

for iSim = 1:nBurst
  totOFDMMuls = 0; 
  totProcMuls = 0;

  for iBWP = BWPs 
    % Input data on active subcarries. Shift around DC
    Xact{iBWP} = 1/sqrt(2)*...
      (sign(randn(nAct(iBWP), bOFDM(iBWP)))+1j*sign(randn(nAct(iBWP), bOFDM(iBWP))));
  
  end

  switch TXproc
   case {'FC-F-OFDM'}
    %% FC-F-OFDM TX PROCESSING
    %% OFDM TX Processing
    for iBWP = BWPs 
      % The data is shifted around zero frequency and FC processing modulates 
      %   the waveform to desired frequency 
      X{iBWP} = shiftData(Xact{iBWP}, L_OFDM_LR(iBWP), nAct(iBWP), bOFDM(iBWP));

      % OFDM modulation and CP insertion. Here, CP-OFDM waveform is created 
      %  at low rate and FC processing interpolates the waveform at high output rate
      x_cp{iBWP} = sqrt(L_OFDM_LR(iBWP))* ...
          OFDMmod(X{iBWP}, nAct(iBWP), bOFDM(iBWP), L_OFDM_LR(iBWP), L_CP_LR(iBWP,:));
      
      % Zero-padding needed for fast convolution
      x_cp{iBWP} = [zeros(prePad, 1); x_cp{iBWP}; zeros(postPad, 1)]; 
    
      % Cumulative complexity for OFDM processing
      totOFDMMuls = totOFDMMuls + bOFDM(iBWP)*FFTcomp(L_OFDM_LR(iBWP));
    end
    
    %% Fast-Convolution TX Processing (synthesis filterbank)
    w_buff = zeros(N, nBlk);
    for iBWP = BWPs     
      % Buffering into overlapping blocks 
      x_buff{iBWP} = buffer(x_cp{iBWP}(:), L, L-Ls, 'nodelay'); 
      
      % Phase rotation in order to maintain the phase continuation between processing 
      %   blocks in the case of frequency-domain shift by Bin0 bins (modulation)
      x_buff{iBWP} = bsxfun(@times, x_buff{iBWP}(:, 1:nBlk), ...
                            exp(2j*pi*Bin0(iBWP)*(0:nBlk-1)*Ls/L));
      
      % FC processing: (see, e.g, [2] for details)
      %   FFT -> FFT-shift -> windowing -> modulation -> IFFT -> discard overlap
      % The high-rate waveform containing both BWPs is combined in frequency domain 
      w_buff = w_buff + M_LN{iBWP}*win{iBWP}*fftshift(fft(x_buff{iBWP}, [], 1), 1);
    
      % Frequency-domain windowing (2*2*nTB multiplications)
      totProcMuls = totProcMuls + nBlk*4*nTB;

      % Cumulative complexity for FC processing (Short FFTs of length L)
      totProcMuls = totProcMuls + nBlk*FFTcomp(L);

    end

    % Conversion to time domain
    buff = sqrt(N/L)*ifft(w_buff);

    % Cumulative complexity for FC processing (Long IFFTs of length N)
    totProcMuls = totProcMuls + nBlk*FFTcomp(N);
    
    % Discard overlap and concatenate the processed blocks
    y = [zeros(Nl,1); vec(S_Ns*buff); zeros(Nl,1)];
      
   case {'WOLA'}
    %% WOLA TX PROCESSING
    y = zeros(N/L*nTotSamp, 1);

    % Here we are using longer extension than possible with EVM requirements 
    % given in TS38.104
    WOLAext = L_CP_HR([1 2], 2)/2; 
    for iBWP = BWPs     
      % Input data on active subcarries. Shift around DC
      X{iBWP} = shiftData(Xact{iBWP}, L_OFDM_HR(iBWP), nAct(iBWP), bOFDM(iBWP));

      % OFDM modulation and CP insertion
      x_cp{iBWP} = sqrt(L_OFDM_HR(iBWP))* ...
          WOLAmod(X{iBWP}, nAct(iBWP), bOFDM(iBWP), L_OFDM_HR(iBWP), ...
                  L_CP_HR(iBWP, :), WOLAext(iBWP));

      % Cumulative complexity for OFDM processing
      totOFDMMuls = totOFDMMuls + bOFDM(iBWP)*FFTcomp(L_OFDM_HR(iBWP));

      % Cumulative complexity for WOLA TX (2*2*2*WOLAext multiplications per symbol)
      totProcMuls = totProcMuls + bOFDM(iBWP)*2*2*2*WOLAext(iBWP);

      % Zero-padding (in order to make the delays of both processing alternatives 
      %   the same) 
      w = [zeros(N/L*prePad-WOLAext(iBWP), 1); ...
           x_cp{iBWP}; ...
           zeros(N/L*postPad-WOLAext(iBWP), 1)]; 

      % Modulation to desired center frequency
      y = y + w.*exp(j*2*pi*Bin0(iBWP).*BS_FC/Fs*(0:N/L*nTotSamp-1).');
   end
    
   case {'CP-OFDM'}
    %% BASIC OFDM TX PROCESSING
    y = zeros(N/L*nTotSamp, 1);

    for iBWP = BWPs     
      % Input data on active subcarries. Shift around DC
      X{iBWP} = shiftData(Xact{iBWP}, L_OFDM_HR(iBWP), nAct(iBWP), bOFDM(iBWP));

      % OFDM modulation and CP insertion
      x_cp{iBWP} = sqrt(L_OFDM_HR(iBWP))* ...
         OFDMmod(X{iBWP}, nAct(iBWP), bOFDM(iBWP), L_OFDM_HR(iBWP), L_CP_HR(iBWP, :));

      % Cumulative complexity for OFDM processing
      totOFDMMuls = totOFDMMuls + bOFDM(iBWP)*FFTcomp(L_OFDM_HR(iBWP));

      % Zero-padding (in order to make the delays of both processing alternatives 
      %   the same) 
      w = [zeros(N/L*prePad, 1); x_cp{iBWP}; zeros(N/L*postPad, 1)]; 

      % Modulation to desired center frequency
      y = y + w.*exp(j*2*pi*Bin0(iBWP).*BS_FC/Fs*(0:N/L*nTotSamp-1).');
   end
    
   case {'TD-F-OFDM'}
    %% TIME-DOMAIN FILTERED OFDM TX PROCESSING
    y = zeros(N/L*nTotSamp, 1);

    for iBWP = BWPs     
      % Input data on active subcarries. Shift around DC
      X{iBWP} = shiftData(Xact{iBWP}, L_OFDM_HR(iBWP), nAct(iBWP), bOFDM(iBWP));

      % OFDM modulation and CP insertion
      x_cp{iBWP} = sqrt(L_OFDM_HR(iBWP))* ...
         OFDMmod(X{iBWP}, nAct(iBWP), bOFDM(iBWP), L_OFDM_HR(iBWP), L_CP_HR(iBWP, :));

      % Cumulative complexity for OFDM processing
      totOFDMMuls = totOFDMMuls + bOFDM(iBWP)*FFTcomp(L_OFDM_HR(iBWP));

      % Zero-padding (to make the delays of all processing alternatives the same) 
      w = [zeros(N/L*prePad, 1); x_cp{iBWP}; zeros(N/L*postPad, 1)]; 

      % Filter design, filtering, and delay compensation
      h = hHann(L_OFDM_HR(iBWP), nAct(iBWP), 14/SCS_fac(iBWP));
      w = conv(h, w);
      
      % Cumulative complexity for filtering processing:
      % Here, real-valued filter is used, therefore, the coefficient symmetry
      % can be utilized. Factor of two due to the real-to-complex multiplies
      filtCmpx = 2*numel(h)/2*numel(w);
      totProcMuls = totProcMuls + 2*numel(h)/2*numel(w); % 
      
      w = w((numel(h)-1)/2+(1:N/L*nTotSamp));
      
      % Modulation to desired center frequency
      y = y + w.*exp(j*2*pi*Bin0(iBWP).*BS_FC/Fs*(0:N/L*nTotSamp-1).');
    
      % Cumulative complexity for modulation:
      % Complex output is multiplied by the complex modulating
      % sequence. Here, four real multiplies are used for complex-to-complex
      % multiplications
      moduCmpx = 4*(N/L*nSamp + numel(h)-1);
      totProcMuls = totProcMuls + moduCmpx;
      
    end % END for iBWP     
    
  end
  
  % Power-spectral density estimate
  [Ysq, freq] = PSDeval(y, Fs, N/L*nSamp, nGridPoints); 
  PSD = PSD + Ysq/nBurst; 
  
  %% TIME-DOMAIN FILTERING RX PROCESSING  
  for iBWP = BWPs
    % Filter design 
    h = hHann(L_OFDM_HR(iBWP), nAct(iBWP), 14/SCS_fac(iBWP));
    delay = (numel(h)-1)/2;

    % Demodulation back to zero frequency 
    z = y(1:N/L*nSamp+pre0+Nl+delay).* ...
        exp(-j*2*pi*Bin0(iBWP).*BS_FC/Fs*((0:N/L*nSamp-1+pre0+Nl+delay)).');
    
    % Filtering
    z = conv(h, z);
      
    % Delay compensation
    z = z(delay+(1:N/L*nSamp+pre0+Nl));
      
    % Remove zero-padding
    z = z(pre0+Nl+(1:N/L*nSamp));  

    % OFDM demodulation
    Yfofdm{iBWP} = 1/sqrt(L_OFDM_HR(iBWP))* ...
        OFDMdem(z, bOFDM(iBWP), L_OFDM_HR(iBWP), L_CP_HR(iBWP, :));
      
    % MSE error on active subcarriers
    MSE_fofdm{iBWP} = MSE_fofdm{iBWP} + ...
        MSE(Yfofdm{iBWP}(iActOFDM{iBWP}, :), X{iBWP}(iActTX{iBWP}, :), bOFDM(iBWP));  
  
  end
  
  %% BASIC OFDM RX PROCESSING  
  for iBWP = BWPs
    % Remove zero-padding
    z = y(pre0+Nl+(1:N/L*nSamp));  

    % Demodulation back to zero frequency 
    z = z.*exp(-j*2*pi*Bin0(iBWP).*BS_FC/Fs*((0:N/L*nSamp-1)+pre0+Nl).');

    % OFDM demodulation
    Yofdm{iBWP} = 1/sqrt(L_OFDM_HR(iBWP))* ...
        OFDMdem(z, bOFDM(iBWP), L_OFDM_HR(iBWP), L_CP_HR(iBWP, :));
      
    % MSE error on active subcarriers
    MSE_ofdm{iBWP} = MSE_ofdm{iBWP} + ...
        MSE(Yofdm{iBWP}(iActOFDM{iBWP}, :), X{iBWP}(iActTX{iBWP}, :), bOFDM(iBWP));  
  
  end
  
  %% WOLA RX PROCESSING  
  for iBWP = BWPs
    % Remove zero-padding
    z = y(pre0+Nl+(1:N/L*nSamp));  

    % Demodulation back to zero frequency 
    z = z.*exp(-j*2*pi*Bin0(iBWP).*BS_FC/Fs*((0:N/L*nSamp-1)+pre0+Nl).');

    % OFDM demodulation
    Ywola{iBWP} = 1/sqrt(L_OFDM_HR(iBWP))* ...
        WOLAdem(z, bOFDM(iBWP), L_OFDM_HR(iBWP), L_CP_HR(iBWP, :));
      
    % MSE error on active subcarriers
    MSE_wola{iBWP} = MSE_wola{iBWP} + ...
        MSE(Ywola{iBWP}(iActOFDM{iBWP}, :), X{iBWP}(iActTX{iBWP}, :), bOFDM(iBWP));  
  
  end
  
  %% FAST-CONVOLUTION RX PROCESSING  
  for iBWP = BWPs
    %% Fast-Convolution RX Processing (analysis filterbank)
    % Zero-padding of input data and buffering to overlapping processing blocks    
    y_buff = buffer(y, N, N-Ns, 'nodelay');

    % FC processing: (see, e.g, [2] for details)
    % FFT-> demodulation -> windowing-> IFFT-shift-> IFFT-> discard overlap
    % win_rx{iBWP} = eye(L);
    z_buff = S_Ls*ifft(ifftshift(win{iBWP}*M_LN{iBWP}.'*fft(y_buff,[], 1), 1), [], 1);
   
    % Phase rotation in order to maintain the phase continuation between processing
    %   blocks in the case of frequency-domain shift by Bin0 bins (modulation)
    z_buff = bsxfun(@times,z_buff(:, 1:nBlk), exp(-2j*pi*(0:nBlk-1)*Bin0(iBWP)*Ls/L));

    % Concatenation of the resulting processing blocks (parallel to serial)
    out = z_buff(:);
    
    %% OFDM RX Processing
    Yfc{iBWP} = 1/sqrt(L_OFDM_HR(iBWP))*...
     OFDMdem(out(pre0*L/N+(1:nSamp)), bOFDM(iBWP), L_OFDM_LR(iBWP), L_CP_LR(iBWP, :));

    % MSE error on active subcarriers
    MSE_fc{iBWP} = MSE_fc{iBWP} + ...
        MSE(Yfc{iBWP}(iActFC{iBWP}, :), X{iBWP}(iActTX{iBWP}, :), bOFDM(iBWP));
    
  end % end for iBWP
  
end % end for nBurst

% Out-of-band emissions
iStopband = find(freq <= -5e6 | freq >= 5e6);
minRipple = max(abs(PSD(iStopband)));

% Scale by number of bursts 
MSE_ofdm{1} = MSE_ofdm{1}/nBurst; MSE_ofdm{2} = MSE_ofdm{2}/nBurst; 
MSE_wola{1} = MSE_wola{1}/nBurst; MSE_wola{2} = MSE_wola{2}/nBurst; 
MSE_fofdm{1} = MSE_fofdm{1}/nBurst; MSE_fofdm{2} = MSE_fofdm{2}/nBurst; 
MSE_fc{1} = MSE_fc{1}/nBurst; MSE_fc{2} = MSE_fc{2}/nBurst; 

Att = -10*log10(minRipple);
MSE_ofdm_bwp1 = 10*log10(mean(MSE_ofdm{1}));
MSE_ofdm_bwp2 = 10*log10(mean(MSE_ofdm{2}));
MSE_fofdm_bwp1 = 10*log10(mean(MSE_fofdm{1}));
MSE_fofdm_bwp2 = 10*log10(mean(MSE_fofdm{2}));
MSE_wola_bwp1 = 10*log10(mean(MSE_wola{1}));
MSE_wola_bwp2 = 10*log10(mean(MSE_wola{2}));
MSE_fc_bwp1   = 10*log10(mean(MSE_fc{1}));
MSE_fc_bwp2   = 10*log10(mean(MSE_fc{2}));
disp([TXproc ':     attenuation at channel edge: ' num2str(Att) ' dB'])
disp([TXproc ':     complexity: OFDM processing: ' int2str(totOFDMMuls) ...
      ' real multiplications'])
disp([TXproc ': complexity: waveform processing: ' int2str(totProcMuls) ...
      ' real multiplications'])
disp([TXproc ':                total complexity: ' ...
      int2str(totOFDMMuls+totProcMuls) ' real multiplications'])
disp(['  CP-OFDM RX: Mean MSE on BPW #1: ' sprintf('%4.1f', MSE_ofdm_bwp1) ' dB'])
disp(['  CP-OFDM RX: Mean MSE on BPW #2: ' sprintf('%4.1f', MSE_ofdm_bwp2) ' dB'])
disp(['     WOLA RX: Mean MSE on BPW #1: ' sprintf('%4.1f', MSE_wola_bwp1) ' dB'])
disp(['     WOLA RX: Mean MSE on BPW #2: ' sprintf('%4.1f', MSE_wola_bwp2) ' dB'])
disp(['FC-F-OFDM RX: Mean MSE on BPW #1: ' sprintf('%4.1f', MSE_fc_bwp1) ' dB'])
disp(['FC-F-OFDM RX: Mean MSE on BPW #2: ' sprintf('%4.1f', MSE_fc_bwp2) ' dB'])
disp(['   F-OFDM RX: Mean MSE on BPW #1: ' sprintf('%4.1f', MSE_fofdm_bwp1) ' dB'])
disp(['   F-OFDM RX: Mean MSE on BPW #2: ' sprintf('%4.1f', MSE_fofdm_bwp2) ' dB'])

%% PLOTTING
figure(1)

if ~isOctave
  colorMap =  get(groot, 'defaultAxesColorOrder');

else
  colorMap =  get(gcf, 'defaultAxesColorOrder');

end

i2s = @(x) int2str(x); % Shorthand for int2str 

% Plot frequency-domain windows of the bandwidth parts
subplot(212)
hline = plot(-L/2:L/2-1, diag(win{1}), 'o-', -L/2:L/2-1, diag(win{2}), 'o-');
set(hline(1), 'MarkerFaceColor', colorMap(1,:))
set(hline(2), 'MarkerFaceColor', colorMap(2,:))
axis([-L/2 L/2-1 -0.1 1.1]), grid on
title(['Frequency-domain window of size L = ' int2str(L)])
xlabel('Subcarrier n'), ylabel('Magnitude')

if strcmpi(TXproc, 'FC-F-OFDM')
  text(-L/2+2, 0.85, ['OFDM transf. length L\_OFDM\_LR = [' i2s(L_OFDM_LR.') ']'])
  text(-L/2+2, 0.75, ['FC-processing short transf. L = ' i2s(L)])
else
  text(-L/2+2, 0.85, ['OFDM transf. length L\_OFDM\_HR = [' i2s(L_OFDM_HR.') ']'])
  text(-L/2+2, 0.75, ['FC-processing short transf. L = ' i2s(L)])
end
text(-L/2+2, 0.65, ['FC-processing long transf. N = ' i2s(N)])
text(-L/2+2, 0.55, ...
     ['OFDM subcarrier spacing SCS\_OFDM = [' i2s(SCS_OFDM.'/1e3) '] kHz'])
text(-L/2+2, 0.45, ['FC-proc. bin spacing BS\_FC = ' i2s(BS_FC.'/1000) ' kHz'])
bwp1s = ['BWP 1: ' i2s(nAct(1)/12) ' PRBs with ' i2s(SCS_OFDM(1)/1e3) ' kHz SCS'];
bwp2s = ['BWP 2: ' i2s(nAct(2)/12) ' PRBs with ' i2s(SCS_OFDM(2)/1e3) ' kHz SCS'];

% Power spectral density of the TX waveform
subplot(211)
yLim = [-100 20]; 
plot([-5 -5 5 5], [yLim+[-1 1] yLim(2:-1:1)+[1 -1]], 'LineWidth', 2), hold on
hfill(1) = fill(nAct(1)/2*[-1 -1 1 1]*SCS_OFDM(1)/1e6+Bin0(1)*BS_FC/1e6, ...
                [yLim+[-1 1] yLim(2:-1:1)+[1 -1]], colorMap(1,:));
hfill(2) = fill([nAct(2)/2*[-1 -1 1 1]]*SCS_OFDM(2)/1e6+Bin0(2)*BS_FC/1e6, ...
                [yLim+[-1 1] yLim(2:-1:1)+[1 -1]], colorMap(2,:));

set(hfill, 'FaceAlpha', 0.3, 'EdgeColor', 'none')
w1 = N/L*(((0:numel(PSD)-1)/numel(PSD))-0.5)*L_OFDM_LR(1);
h = plot(w1*SCS_OFDM(1)/1e6, 10*log10(PSD));
grid on, hold off
title(['PSD of the ' TXproc ' signal, f0 = [' ...
       sprintf('%5.2f', Bin0*BS_FC/1e6) '] MHz'])
ylim(yLim), grid on, xlim(Fs/2*[-1 1]/1e6),
legend('10 MHz channel', bwp1s, bwp2s, [TXproc ' waveform'], 'Location', 'South')
xlabel('Frequency in MHz'), ylabel('PSD in dB')
set(gca, 'GridAlpha', 1, 'GridColor', 0.85*[1 1 1])

%% Mean-squared errors on bandwidth parts
% Basic OFDM
figure(2), axf2 = subplot(121);
hline = plot(-nAct(1)/2:nAct(1)/2-1, 10*log10(MSE_ofdm{1}), 'o-', ...
             -nAct(2)/2:nAct(2)/2-1, 10*log10(MSE_ofdm{2}), 'o-');
set(hline, 'MarkerFaceColor', 'w')
ylabel('MSE error in dB'); xlabel('Active subcarrier n'), grid on
title('Basic OFDM RX: MSE error')
yLim1 = ylim;
legend(bwp1s, bwp2s, 'Location', 'South')
set(gca, 'GridAlpha', 1, 'GridColor', 0.15*[1 1 1])

% WOLA
figure(3), axf3 = subplot(121);
hline = plot(-nAct(1)/2:nAct(1)/2-1, 10*log10(MSE_wola{1}), 'o-', ...
             -nAct(2)/2:nAct(2)/2-1, 10*log10(MSE_wola{2}), 'o-');
set(hline, 'MarkerFaceColor', 'w')
ylabel('MSE error in dB'); xlabel('Active subcarrier n'), grid on
title('WOLA RX: MSE error')
yLim2 = ylim;
legend(bwp1s, bwp2s, 'Location', 'North')
set(gca, 'GridAlpha', 1, 'GridColor', 0.15*[1 1 1])

% FC-F-OFDM
figure(4), axf4 = subplot(121);
hline = plot(-nAct(1)/2:nAct(1)/2-1, 10*log10(MSE_fc{1}), 'o-', ...
             -nAct(2)/2:nAct(2)/2-1, 10*log10(MSE_fc{2}), 'o-');
set(hline, 'MarkerFaceColor', 'w')
ylabel('MSE error in dB'); xlabel('Active subcarrier n'), grid on
title('FC-F-OFDM RX: MSE error')
yLim3 = ylim;
legend(bwp1s, bwp2s, 'Location', 'North')
set(gca, 'GridAlpha', 1, 'GridColor', 0.15*[1 1 1])

% TD-F-OFDM
figure(5), axf5 = subplot(121);
hline = plot(-nAct(1)/2:nAct(1)/2-1, 10*log10(MSE_fofdm{1}), 'o-', ...
             -nAct(2)/2:nAct(2)/2-1, 10*log10(MSE_fofdm{2}), 'o-');
set(hline, 'MarkerFaceColor', 'w')
ylabel('MSE error in dB'); xlabel('Active subcarrier n'), grid on
title('Time-domain filtered OFDM RX: MSE error')
yLim3 = ylim;
minYlim = min([yLim1(1) yLim2(1) yLim3(1)]);
maxYlim = max([yLim1(2) yLim2(2) yLim3(2)]);
set([axf2 axf3 axf4 axf5], 'ylim', [floor(minYlim/5)*5 ceil(maxYlim/5)*5])
legend(bwp1s, bwp2s, 'Location', 'North')
set(gca, 'GridAlpha', 1, 'GridColor', 0.15*[1 1 1])
  
iBWP = 2;
if numel(BWPs) == 1
  iBWP = 1;
end

%% Constellations on bandwidth part # 1
% Basic OFDM
figure(2), subplot(122)
r = plot(Yofdm{iBWP}(iActOFDM{iBWP}, :), '+', 'Color', colorMap(5,:)); hold on
t = plot(X{iBWP}(iActTX{iBWP}, :), 'x', 'Color', colorMap(1,:)); hold off, 
axis('equal'), grid on
title('Basic OFDM RX: Constellation'); 
legend([r(1) t(1)], 'RX{ }', 'TX{ }', 'Location', 'East')
xlabel('Real part'); ylabel('Imaginary part'); 
set(gca, 'GridAlpha', 1, 'GridColor', 0.15*[1 1 1])

% WOLA
figure(3), subplot(122)
r = plot(Ywola{iBWP}(iActOFDM{iBWP}, :), '+', 'Color', colorMap(5,:)); hold on
t = plot(X{iBWP}(iActTX{iBWP}, :), 'x', 'Color', colorMap(1,:)); hold off, 
axis('equal'), grid on
title('WOLA RX: Constellation'); 
legend([r(1) t(1)], 'RX{ }', 'TX{ }', 'Location', 'East')
xlabel('Real part'); ylabel('Imaginary part'); 
set(gca, 'GridAlpha', 1, 'GridColor', 0.15*[1 1 1])

% FC-F-OFDM
figure(4), subplot(122)
r = plot(Yfc{iBWP}(iActFC{iBWP}, :), '+', 'Color', colorMap(5,:)); hold on
t = plot(X{iBWP}(iActTX{iBWP}, :), 'x', 'Color',  colorMap(1,:)); 
hold off, axis('equal'), grid on
title('FC-F-OFDM RX: Constellation'); 
legend([r(1) t(1)], 'RX{ }', 'TX{ }', 'Location', 'East')
xlabel('Real part'); ylabel('Imaginary part'); 
set(gca, 'GridAlpha', 1, 'GridColor', 0.15*[1 1 1])

% TD-F-OFDM
figure(5), subplot(122)
r = plot(Yfofdm{iBWP}(iActOFDM{iBWP}, :), '+', 'Color', colorMap(5,:)); hold on
t = plot(X{iBWP}(iActTX{iBWP}, :), 'x', 'Color', colorMap(1,:)); hold off, 
axis('equal'), grid on
title('Time-domain filtered RX: Constellation'); 
legend([r(1) t(1)], 'RX{ }', 'TX{ }', 'Location', 'East')
xlabel('Real part'); ylabel('Imaginary part'); 
set(gca, 'GridAlpha', 1, 'GridColor', 0.15*[1 1 1])

end

%% SUBFUNCTIONS
% ------------------------------------------------------------------------------------
function mse = MSE(Y, X, bOFDM)

  % Mean-squared error over active subcarriers

  mse = sum(abs(Y-X).^2, 2)/bOFDM;

end

% ------------------------------------------------------------------------------------
function X = shiftData(Xact, L_OFDM, nAct, bOFDM)

  % Shift data around zero frequency. 
  
  X = circshift([fftshift(Xact, 1); zeros(L_OFDM-nAct, bOFDM)], [-floor(nAct/2) 0]); 
  
end  
  
% ------------------------------------------------------------------------------------
function X = mapData(Xact, L_OFDM, nAct, bOFDM, shift)

  % Shift data around zero frequency. 
  
  X = circshift([Xact; zeros(L_OFDM-nAct, bOFDM)], [shift 0]); 
  
end  
  
% ------------------------------------------------------------------------------------
function x_cp = OFDMmod(X, nAct, bOFDM, L_OFDM, L_CP)

  % OFDM modulation and the inclucion of CP

  x_cp = zeros(bOFDM*L_OFDM+sum(L_CP), 1);
  for iSymb = 1:bOFDM
    x = ifft(X(:, iSymb));
    
    x_cp((iSymb-1)*L_OFDM + sum(L_CP(1:iSymb-1)) + (1:L_OFDM+L_CP(iSymb))) = ...
        [x(L_OFDM-L_CP(iSymb)+1:L_OFDM, :); x]; % Add cyclic prefix
    
  end

end

% ------------------------------------------------------------------------------------
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
    time_offset = L_CP(iSymb)/2; 
  
    OFDMsymb = z(L_OFDM*(iSymb-1)+sum(L_CP(1:iSymb))-time_offset+(1:L_OFDM));
    
    Y(:, iSymb) = fft(circshift(OFDMsymb, [-time_offset 0]));    
  
  end     

end

% ------------------------------------------------------------------------------------
function x_cp = WOLAmod(X, nAct, bOFDM, L_OFDM, L_CP, WOLAext)

  % WOLA TX processing

  % Function for generating the transition band weights
  tapFun = @(T, root)((cos(2*pi*[-(T-0):-1]/(2*(T+1)))+1)/2).'.^(1/root);
 
  tap = tapFun(2*WOLAext, 1); % Window taper
  
  x_cp = zeros(bOFDM*L_OFDM+sum(L_CP) + 2*WOLAext, 1);    
  for iSymb = 1:bOFDM
    WOLAwin = [tap; ones(L_OFDM+L_CP(iSymb)-2*WOLAext, 1); tap(2*WOLAext:-1:1)];
    
    OFDMsymb = ifft(X(:, iSymb));

    % Add postfix and prefix
    WOLAsymb = [
        OFDMsymb(L_OFDM-L_CP(iSymb)-WOLAext+1:L_OFDM, :); % CP + EXT samples from end
        OFDMsymb; 
        OFDMsymb(1:WOLAext) % EXT samples from the begin
               ];
    
    windowedSymb = WOLAsymb.*WOLAwin;
    
    inx = (iSymb-1)*L_OFDM + sum(L_CP(1:iSymb-1)) + (1:L_OFDM+L_CP(iSymb)+2*WOLAext);
    x_cp(inx) = x_cp(inx) + windowedSymb;

  end
  
  C_WOLA = 2*2*numel(tap)*bOFDM;

end

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

  % WOLA RX processing

  WOLAext = L_CP(2)/2; 
  
  % Function for generating the transition band weights
  tapFun = @(T, root)((cos(2*pi*[-(T-0):-1]/(2*(T+1)))+1)/2).'.^(1/root);
 
  tap = tapFun(2*WOLAext, 1);
  
  WOLAwin = [tap; ones(L_OFDM-2*WOLAext, 1); tap(2*WOLAext:-1:1)];
    
  Y = zeros(L_OFDM, bOFDM); 
  for iSymb = 1:bOFDM
    iCPOFDM = L_OFDM*(iSymb-1)+sum(L_CP(1:iSymb-1))+L_CP(iSymb)+(1-2*WOLAext:L_OFDM);
    
    % Windowing 
    windowedSymb = WOLAwin.*z(iCPOFDM);
    
    % Indices of the overlapping parts
    inx1 = WOLAext + (1:WOLAext);
    inx2 = 2*WOLAext + L_OFDM-WOLAext + (1:WOLAext);
    inx3 = 2*WOLAext + (1:L_OFDM-2*WOLAext);
    inx4 = (1:WOLAext);
    inx5 = WOLAext+L_OFDM-WOLAext+(1:WOLAext);
    
    % Overlap-and-add 
    WOLAsymb = [
        windowedSymb(inx1) + windowedSymb(inx2); ...
        windowedSymb(inx3);
        windowedSymb(inx4) + windowedSymb(inx5); ...    
        ];
    
    % Compensate for delay and demodulate
    Y(:, iSymb) = fft(circshift(WOLAsymb, [-WOLAext 0]));    
  
  end     

end

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

  % Frequency-domain window function
  
  % Optimized transition-band weights. Here, the optimization 
  % is carried out following the techniques given in [2]
  tap = [ % JYK; 27.11.19, 14:25
          -0.000340952946487
          +0.000769433263641
          +0.003043703626772
          -0.012311509605436
          -0.052593175378407
          -0.081557396098516
          -0.018092949343188
          +0.183410321328965
          +0.481282119474811
          +0.759888460602338
          +0.929330128417243
          +0.991705584048036
        ];

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

  win = sparse(1:L, 1:L, 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 h = hHann(L_OFDM, nAct, excessBW)

  % Filter design for time-domain filtered OFDM as proposed in [1]. 
  % [1] X. Zhang, M. Jia, L. Chen, J. Ma, and J. Qiu, "Filtered-OFDM - Enabler for 
  % flexible waveform in the 5th generation cellular networks," in IEEE Global 
  % Communications Conference (GLOBECOM), Dec. 2015, pp. 16.

  L = floor(L_OFDM/2);    % Filter order
  n = (0:L)-floor(L/2);
  
  h = sinc(n*(nAct+excessBW)/L_OFDM);   % Rectangular in frequency -> sinc in time  
  w = (0.5*(1+cos(2*pi*n/(L-1)))).^0.5; % Hann window
  
  h = h.*w;     % Windowing
  h = h/sum(h); % Normalization

end
% ------------------------------------------------------------------------------------
function multSR = FFTcomp(N)

  % Complexity of split-radix for power-of-two length N
  log2N = log2(N);
  if abs(mod(log2N, 1)) > 100*eps
    error('This complexity applies only for power-of-two transform sizes')
  
  end
  multSR = N*(log2N - 3) + 4;

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

Simulation Results

>> filteredOFDM('FC-F-OFDM')
FC-F-OFDM:     attenuation at channel edge: 85.7839 dB
FC-F-OFDM:     complexity: OFDM processing: 79016 real multiplications
FC-F-OFDM: complexity: waveform processing: 429440 real multiplications
FC-F-OFDM:                total complexity: 508456 real multiplications
  CP-OFDM RX: Mean MSE on BPW #1: -33.0 dB
  CP-OFDM RX: Mean MSE on BPW #2: -29.9 dB
     WOLA RX: Mean MSE on BPW #1: -54.5 dB
     WOLA RX: Mean MSE on BPW #2: -42.8 dB
FC-F-OFDM RX: Mean MSE on BPW #1: -51.8 dB
FC-F-OFDM RX: Mean MSE on BPW #2: -45.2 dB
   F-OFDM RX: Mean MSE on BPW #1: -51.6 dB
   F-OFDM RX: Mean MSE on BPW #2: -45.1 dB