%% Simplified model for the generalized 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 some non-trivial weights
tap = [ 0.4463042 0.2183986 0.6073331 0.8429762 1.2444025 1.6888477 1.7148595 0.9621689]; Ntap = numel(tap); 
a_m = [-0.3306919 1.3273771 1.6471274 1.1502148 0.7189577 0.7133351 0.9474267  1.1462551 ...
	1.2006866 1.1075530 0.8985521 0.7163255 0.8006488 1.2200905 1.5733476  1.1707125 ];
A = sparse(diag(repmat(a_m(:), [L/numel(a_m), 1]))); S = sparse(diag(hanning(N, 'periodic'))); 
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 % Loop for bursts
  %% TX PROCESSING  ========================================
  w = zeros((bOFDM(1)+1)*N, 1); % Pre-allocate output vector
  for k = 1:2 % Loop for subbands    
    inx = (-L/L_OFDM(k)*(nAct(k)/2+1)-Ntap:L/L_OFDM(k)*(nAct(k)/2+1)-1+Ntap)+L/2;       % Optimized 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_ZP = [zeros(L/2, 1); vec(sqrt(L_OFDM(k))*ifft(X{k})); zeros(L/2, 1)]; % Eq. (7c)    

    if (1 == 1) % Processing model with block diagonal synthesis matrix 
      % Here Fm is generated as a sum of two block diagonal matrices without overlap
      tmp = cell(ceil(Rm/2), 1); [tmp{:}] = deal(sparse(Fmr)); 
      F = blkdiag(tmp{:}) + blkdiag(Fmr(N/2+(1:N/2),L/2+(1:L/2)),tmp{1:floor(Rm/2)},Fmr(1:N/2,1:L/2)); % Eq. (7b)
      % 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_ZP;
    else % Alternative processing model with overlapping input blocks
      x_buff = buffer(x_ZP, L, L/2, 'nodelay'); % Buffering of input data Eq. (15a)
      for r = 0:Rm-1                            % Overlap-and-add processing
	w(r*N/2+(1:N)) = w(r*N/2+(1:N)) + sqrt(N/L)*Fmr*x_buff(:, r+1);
      end
    end

  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')