% Simplified model for the basic FC processing, (c) Juha Yli-Kaakinen, 6.9.2019
clear all, ran_seed = 10; % rng(ran_seed);
vec = @(x) x(:);   % Vectorization function 
L = 128;           % Short FC processing transform size
N = 1024;          % Long FC processing transform size
Fs = 15e3*N;       % Sampling frequency
Bin0 = [448 576];  % Center bins of the subbands
L_OFDM = [128 64]; % OFDM transform sizes (15 kHz and 30 kHz SCSs for 15.36 MHz sampling rate)
bOFDM  = [14 28];  % Number of OFDM symbols for subbands (14 for 15 kHz SCS and 28 for 30 kHz SCS)
nAct = 12*[2 1];   % Number of active subcarriers on each subband
W_L = fft(eye(L)); WI_L = ifft(eye(L)); WI_N = ifft(eye(N));               % FFT and IFFT matrices
P = sparse(fftshift(eye(L), 1));                                           % FFT shift
M_LN = @(c) sparse(mod(ceil(-L/2:L/2-1)+c, N)+1, 1:L, ones(1, L), N, L);   % Frequency-domain modulation

% Windowing matrices with analytical weights
A = sparse(diag(ones(L, 1))); S = sparse(diag([zeros(N/4, 1); ones(N/2, 1); zeros(N/4, 1)])); % Eq. (15)
Ntap = 8; tap = sqrt(hanning(2*Ntap+1, 'symmetric')); tap = tap(1:Ntap).'; % RRC transition-band bins (8 taps)

nSim = 100; PSD = zeros((bOFDM(1)+1)*N, 1); MSE{1} = zeros(nAct(1), 1); MSE{2} = zeros(nAct(2), 1); % Preallocate MSE and PSD
Rm = 2*bOFDM(1)*L_OFDM(1)/L(1)+1; % Number of processing blocks for this trivial case
for iSim = 1:nSim 
  
  %% TX PROCESSING  ========================================
  w = zeros((bOFDM(1)+1)*N, 1); % Pre-allocate output vector
  for k = 1:2 % Loop for subbands
    % Analytical subband-wise frequency-domain windows 
    inx = (-L/L_OFDM(k)*(nAct(k)/2+1)-Ntap:L/L_OFDM(k)*(nAct(k)/2+1)-1+Ntap)+L/2;       % Analytical subband-wise 
    D{k} = sparse(inx, inx, [tap ones(1, L/L_OFDM(k)*(nAct(k)+2)) fliplr(tap)], L, L);  %   frequency-domain windows 

    C = WI_L*D{k}*P*W_L;        % Convolution (Eq. 12b)
    U = WI_N*M_LN(Bin0(k))*W_L; % Interpolation (Eq. 12c)
    Fmr = S*U*C*A;              % Overall processing (Eq. 12a)
    
    % Generate random symbols on active subcarriers and allocate around DC bin 
    Xact{k} = 1/sqrt(2)*(sign(randn(nAct(k), bOFDM(k))) + 1j*sign(randn(nAct(k), bOFDM(k))));    
    X{k} = circshift([fftshift(Xact{k}, 1); zeros(L_OFDM(k)-nAct(k), bOFDM(k))], [-floor(nAct(k)/2) 0]);
    
    % OFDM modulation, vectorization, and zero-padding 
    x = [zeros(L/2, 1); vec(sqrt(L_OFDM(k))*ifft(X{k})); zeros(L/2, 1)]; 
    
    % Block-diagonal matrix with overlaping blocks (sum of two block-diagonal matrices without overlap)
    tmp = cell(bOFDM(1), 1); [tmp{:}] = deal(sparse(Fmr)); 
    F = blkdiag(Fmr,tmp{:}) + blkdiag(Fmr(N/2+(1:N/2),L/2+(1:L/2)),tmp{:},Fmr(1:N/2,1:L/2));
    
    % Block-wise processing and combination of the subbands. Combination can also be done in 
    % frequency-domain before the long IFFT and synthesis window as explained in manuscript
    w = w + sqrt(N/L)*F*x;
    
  end
  W = fft(w); PSD = PSD + 1/(N/L*L/2*Rm)*abs(W).^2/nSim; % PSD estimate

  %% RX PROCESSING ========================================
  w_sub = w(N/L*L/2+(1:N/L*bOFDM(1)*L_OFDM(1))); % Remove zero-padding
  
  Y{1} = 1/sqrt(N/L*L_OFDM(1))*fft(reshape(w_sub, N/L*L_OFDM(1), bOFDM(1))); % BASIC OFDM RX PROCESSING for subchannel 1
  iAct_RX = mod([1:nAct(1)]-nAct(1)/2-1+L_OFDM(1)/L*Bin0(1), N/L*L_OFDM(1))+1;
  err{1} = X{1}([L_OFDM(1)-nAct(1)/2+1:L_OFDM(1) 1:nAct(1)/2], 1:bOFDM(1)) - Y{1}([iAct_RX], 1:bOFDM(1));
  MSE{1} = MSE{1} + sum(abs(err{1}).^2, 2)/bOFDM(1)/nSim;
  
  Y{2} = 1/sqrt(N/L*L_OFDM(2))*fft(reshape(w_sub, N/L*L_OFDM(2), bOFDM(2))); % BASIC OFDM RX PROCESSING for subchannel 2
  iAct_RX = mod([1:nAct(2)]-nAct(2)/2-1+L_OFDM(2)/L*Bin0(2), N/L*L_OFDM(2))+1;
  err{2} = X{2}([L_OFDM(2)-nAct(2)/2+1:L_OFDM(2) 1:nAct(2)/2], 1:bOFDM(2)) - Y{2}([iAct_RX], 1:bOFDM(2));
  MSE{2} = MSE{2} + sum(abs(err{2}).^2, 2)/bOFDM(2)/nSim;
  
end

subplot(421), plot(1:L, diag(A)), title('Analysis window'),  xlim([1 L]), xlabel('Samples'), ylabel('Magnitude')
subplot(423), plot(1:N, diag(S)), title('Synthesis window'), xlim([1 N]), xlabel('Samples'), ylabel('Magnitude')
subplot(222), plot(1:L, diag(D{1}), 1:L, diag(D{2}), '--'), 
title('Frequency-domain windows'), xlim([1 L]), xlabel('Samples'), ylabel('Magnitude'), legend('Subband 1', 'Subband 2')
subplot(223), plot(linspace(-1, 1, numel(PSD))*Fs/2/1e6, 10*log10(PSD)), ylim([-120 20]), grid on
title('PSD estimate of the TX waveform'), xlabel('Frequency in MHz'), ylabel('Power in dB'), xlim([-Fs/2 Fs/2]/1e6)
subplot(224), plot(-nAct(1)/2:nAct(1)/2-1, 10*log10(MSE{1}), 'o-', -nAct(2)/2:nAct(2)/2-1, 10*log10(MSE{2}), 'o--', ...
  'MarkerFaceC', 'w'), title('MSEs on active subcarriers'),  xlabel('Active subcarriers'), 
ylabel('MSE in dB'), xlim([-max(nAct)/2 max(nAct)/2-1]), legend('Subband 1', 'Subband 2')