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. DOI: 10.1002/9781119471509.w5GRef002

Abstract

Orthogonal frequency-division multiplexing (OFDM) has been adopted as the waveform of choice in the existing and emerging broadband wireless communication systems for a number of advantages it can offer. Nevertheless, investigations of more advanced multicarrier transmission schemes have continued with the aim of eliminating or mitigating its essential limitations. This article discusses multicarrier schemes with enhanced spectrum localization, which manage to reduce the spectral sidelobes of plain OFDM that are problematic in various advanced communication scenarios. These include schemes for enhancing the OFDM waveform characteristics through additional signal processing as well as filter-bank multicarrier (FBMC) waveforms utilizing frequency-selective filter banks instead of plain (inverse) discrete Fourier transform processing for waveform generation and demodulation.

KEYWORDS

multicarrier; filtered orthogonal frequency-division multiplexing (OFDM); filter bank; filter-bank multicarrier with offset quadrature amplitude modulation (FBMC/OQAM); filtered multitone (FMT); generalized frequency-division multiplexing (GFDM); fast convolution; windowed overlap-and-add (WOLA); universally filtered multicarrier (UFMC); fifth generation (5G) new radio

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 several 5G-NR processing 
%%              alternatives. Supported processing models are plain CP-OFDM, windowed
%%              overlap-and-add (WOLA), time-domain filtered OFDM (TD-F-OFDM), and 
%%              fast-convolution filtered OFDM (FC-F-OFDM) TX and RX processing.
%%              Computational complexity evaluations are given for all models.
%% 
%% 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: Mon Mar 16 18:36:45 2020
%%           By: ylikaaki
%%     Update #: 4695
%% Keywords: 
%% Usage: Call with filteredOFDM('Model') where 'Model' is either 'CP-OFDM, 
%%        'FC-F-OFDM', 'TD-F-OFDM', or 'WOLA'.
%% 
%% 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 (inconjunction 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
%% [3] Juha Yli-Kaakinen, Markku Renfors, and Eleftherios Kofidis, "Filtered 
%%     Multicarrier Transmission," Chapter 2, Wiley 5G Ref, eds. R. Tafazolli, 
%%     C.-L. Wang, and P. Chatzimisios, Wiley, 2020.
%%     https://doi.org/10.1002/9781119471509.w5GRef002
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function filteredOFDM(TXproc)

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

ran_seed = 10; 

isOctave = exist('OCTAVE_VERSION', 'builtin') == 5;
if isOctave 
  pkg load communications
  randn('seed', ran_seed);
  rand('seed', ran_seed);
  colorMap =  get(gcf, 'defaultAxesColorOrder');

else
  rng(ran_seed);
  colorMap =  get(groot, 'defaultAxesColorOrder');

end

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

MOD_ORDER = 4;              % Modulation order 
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/L)
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) + 4;  % 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) + 4; % For the 1st symbols in a slot

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

nTB = 12;                   % Number of transition band weights in FC processing
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-and-save matrices
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;

BWPs = 1:2;         % Bandwidth parts to be processing
nBurst = 100;       % Number of bursts to be transmitted
nGridPoints = 2^16; % Number of grid points used for evaluating the PSD estimates

% Preallocate vectors for MSE and PSD evaluations
for iBWP = BWPs
  disp(['------------ BWP # ' int2str(iBWP) ' ------------'])
  disp([' OFDM length: ' int2str(L_OFDM_HR(iBWP)) ])
  disp(['  CP lengths: [' int2str(L_CP_HR(iBWP, 1:2)) ']'])
  disp(['Burst length: ' int2str(I*nSamp)])

  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
disp(['---------------------------------'])

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;
  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 % END for iBWP

for iSim = 1:nBurst
  % Counters for complexity evaluations
  totOFDMMuls = 0; totProcMuls = 0;

  for iBWP = BWPs 
    % QAM symbols on active subcarries, average energy normalized to unity
    tx_data = randi(MOD_ORDER, nAct(iBWP), bOFDM(iBWP)) - 1;
    Xact{iBWP} = 1/sqrt(2/3*(MOD_ORDER-1))*qammod(tx_data, MOD_ORDER);
  
  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 (OFDM IFFTs)
      totOFDMMuls = totOFDMMuls + bOFDM(iBWP)*FFTcomp(L_OFDM_LR(iBWP));
    
    end % END for iBWP
    
    %% Fast-Convolution TX Processing (synthesis filterbank)
    z_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 
      z_buff = z_buff + M_LN{iBWP}*win{iBWP}*fftshift(fft(x_buff{iBWP}, [], 1), 1);
    
      % Frequency-domain windowing (2*2*nTB multiplications)
      totProcMuls = totProcMuls + 4*nBlk*nTB;

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

    end % END for iBWP

    % Conversion to time domain
    buff = sqrt(N/L)*ifft(z_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
    WOLA2xext = L_CP_HR([1 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, :), WOLA2xext(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 + 2*2*WOLA2xext(iBWP)*bOFDM(iBWP);

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

      % Modulation to desired center frequency
      y = y + w.*exp(j*2*pi*Bin0(iBWP).*BS_FC/Fs*(0:N/L*nTotSamp-1).');
    
    end % END for iBWP
    
   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 % END for iBWP
    
   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 % END switch TXproc
  
  % 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));  

    % Timing offset
    toff = -L_CP_HR(iBWP, 2)/2;
    
    % OFDM demodulation
    Yfofdm{iBWP} = 1/sqrt(L_OFDM_HR(iBWP))* ...
      OFDMdem(z, bOFDM(iBWP), L_OFDM_HR(iBWP), L_CP_HR(iBWP, :), toff);
      
    % MSE error on active subcarriers
    MSE_fofdm{iBWP} = MSE_fofdm{iBWP} + ...
      MSE(Yfofdm{iBWP}(iActOFDM{iBWP}, :), X{iBWP}(iActTX{iBWP}, :), bOFDM(iBWP));  
  
  end % END for iBWP
  
  %% 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).');

    % Timing offset
    toff = -L_CP_HR(iBWP, 2)/2;
    
    % OFDM demodulation
    Yofdm{iBWP} = 1/sqrt(L_OFDM_HR(iBWP))* ...
      OFDMdem(z, bOFDM(iBWP), L_OFDM_HR(iBWP), L_CP_HR(iBWP, :), toff);
      
    % MSE error on active subcarriers
    MSE_ofdm{iBWP} = MSE_ofdm{iBWP} + ...
      MSE(Yofdm{iBWP}(iActOFDM{iBWP}, :), X{iBWP}(iActTX{iBWP}, :), bOFDM(iBWP));  
  
  end % END for iBWP
  
  %% WOLA RX PROCESSING  
  WOLA2xext = L_CP_HR([1 2], 2);     
  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, :), WOLA2xext(iBWP));
      
    % MSE error on active subcarriers
    MSE_wola{iBWP} = MSE_wola{iBWP} + ...
        MSE(Ywola{iBWP}(iActOFDM{iBWP}, :), X{iBWP}(iActTX{iBWP}, :), bOFDM(iBWP));  
  
  end % END for iBWP
  
  %% 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
    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(:);
    
    % Timing offset
    toff = -L_CP_LR(iBWP, 2)/2;
    
    %% 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, :), toff);

    % 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);
maxRipple = max(abs(PSD(iStopband)));

Att = -10*log10(maxRipple);
MSE_ofdm_bwp1  = 10*log10(mean(MSE_ofdm{1}/nBurst));
MSE_ofdm_bwp2  = 10*log10(mean(MSE_ofdm{2}/nBurst));
MSE_fofdm_bwp1 = 10*log10(mean(MSE_fofdm{1}/nBurst));
MSE_fofdm_bwp2 = 10*log10(mean(MSE_fofdm{2}/nBurst));
MSE_wola_bwp1  = 10*log10(mean(MSE_wola{1}/nBurst));
MSE_wola_bwp2  = 10*log10(mean(MSE_wola{2}/nBurst));
MSE_fc_bwp1    = 10*log10(mean(MSE_fc{1}/nBurst));
MSE_fc_bwp2    = 10*log10(mean(MSE_fc{2}/nBurst));
disp([TXproc ' TX:     attenuation at channel edge: ' sprintf('%4.1f', Att) ' dB'])
disp([TXproc ' TX:     complexity: OFDM processing: ' int2str(totOFDMMuls) ...
      ' real multiplications'])
disp([TXproc ' TX: complexity: waveform processing: ' int2str(totProcMuls) ...
      ' real multiplications'])
disp([TXproc ' TX:                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)

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(['FC-processing 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(1), 'FaceColor', [244 203 186]/255, 'EdgeColor', 'none')
set(hfill(2), 'FaceColor', [179 213 235]/255, 'EdgeColor', 'none')
set(gca, 'Layer', 'top')

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 estimate of the ' TXproc ' signal, BWP center frequencies: 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 % END function filteredOFDM

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

  % Mean-squared error over active subcarriers

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

end % END function MSE

% --------------------------------------------------------------------------------------
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 % END function shiftData
  
% --------------------------------------------------------------------------------------
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 % END function mapData
  
% --------------------------------------------------------------------------------------
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 for iSymb

end % END function OFDMmod

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

  % Removal of the CP and OFDM demodulation

  Y = zeros(L_OFDM, bOFDM); 
  for iSymb = 1:bOFDM

    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 for iSymb

end % END function OFDMdem

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

  % 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(WOLA2xext, 1); % Window taper
  
  x_cp = zeros(bOFDM*L_OFDM+sum(L_CP) + WOLA2xext, 1);    
  for iSymb = 1:bOFDM
    WOLAwin = [tap; ones(L_OFDM+L_CP(iSymb)-WOLA2xext, 1); tap(WOLA2xext:-1:1)];
    
    OFDMsymb = ifft(X(:, iSymb));

    % Add postfix and prefix
    WOLAsymb = [
        OFDMsymb(L_OFDM-L_CP(iSymb)-ceil(WOLA2xext/2)+1:L_OFDM, :); 
        OFDMsymb; 
        OFDMsymb(1:floor(WOLA2xext/2))
               ];
    
    windowedSymb = WOLAsymb.*WOLAwin;
    
    inx = (iSymb-1)*L_OFDM + sum(L_CP(1:iSymb-1)) + (1:L_OFDM+L_CP(iSymb)+WOLA2xext);
    x_cp(inx) = x_cp(inx) + windowedSymb;

  end % END for iSymb
  
end % END function WOLAmod

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

  % WOLA RX processing

  % Function for generating the transition band weights
  tapFun = @(T, root)((cos(2*pi*[-(T-0):-1]/(2*(T+1)))+1)/2).'.^(1/root); 
  
  Y = zeros(L_OFDM, bOFDM);
  time_offset = 0;
  
  for iSymb = 1:bOFDM

    tap = tapFun(WOLA2xext, 1);
    
    WOLAwin = [tap; ones(L_OFDM-WOLA2xext, 1); tap(WOLA2xext:-1:1)];
    
    iCPOFDM = L_OFDM*(iSymb-1) + sum(L_CP(1:iSymb-1)) + L_CP(iSymb) + ...
              time_offset + (1-WOLA2xext:L_OFDM);
    
    % Windowing 
    windowedSymb = WOLAwin.*z(iCPOFDM);
    
    % Indices of the overlapping parts
    inx1 = WOLA2xext + (1:L_OFDM-WOLA2xext);
    inx2 = (1:WOLA2xext);
    inx3 = L_OFDM+(1:WOLA2xext);
    
    % Overlap-and-add 
    WOLAsymb = [
        windowedSymb(inx1);
        windowedSymb(inx2) + windowedSymb(inx3); ...    
               ];
    
    % Compensate for delay and demodulate
    Y(:, iSymb) = fft(circshift(WOLAsymb((1:L_OFDM)), [time_offset 0]));  
    
  end     

end % END function WOLAdem

% --------------------------------------------------------------------------------------
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 described 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 % END function genWin

% --------------------------------------------------------------------------------------
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 % END function PSDeval

% --------------------------------------------------------------------------------------
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 % END function hHann

% --------------------------------------------------------------------------------------
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 % END function FFTcomp

% --------------------------------------------------------------------------------------

Simulation Results

octave:1> filteredOFDM('FC-F-OFDM')
------------ BWP # 1 ------------
 OFDM length: 1024
  CP lengths: [80  72]
Burst length: 15360
------------ BWP # 2 ------------
 OFDM length: 512
  CP lengths: [44  36]
Burst length: 15360
---------------------------------
FC-F-OFDM TX:     attenuation at channel edge: 85.7 dB
FC-F-OFDM TX:     complexity: OFDM processing: 79016 real multiplications
FC-F-OFDM TX: complexity: waveform processing: 416020 real multiplications
FC-F-OFDM TX:                total complexity: 495036 real multiplications
  CP-OFDM RX: Mean MSE on BPW 1: -32.8 dB
  CP-OFDM RX: Mean MSE on BPW 2: -29.7 dB
     WOLA RX: Mean MSE on BPW 1: -54.6 dB
     WOLA RX: Mean MSE on BPW 2: -42.7 dB
FC-F-OFDM RX: Mean MSE on BPW 1: -51.8 dB
FC-F-OFDM RX: Mean MSE on BPW 2: -45.0 dB
   F-OFDM RX: Mean MSE on BPW 1: -51.6 dB
   F-OFDM RX: Mean MSE on BPW 2: -44.9 dB