close all;
clear;
clc;
rng('default');

%% ANALYSIS OF OFDM BASED SINGLE CARRIER SYSTEM WITH SINGLE SUBCHANNEL ENABLED WITH Tx FILTERING
% This code calculates the measured SNR, EVM and BER for the single carrier systems with
% filtering at Transmitter.
% We have populated only the middle subchannel in this code.

% The frequency mapping of the single carrier system has a total of 16 subchannels
% with the middle 12 subchannels for data and two subchannels at the extreme end on
% either side for guard.
% In the time domain, for every 9 symbols, the middle symbol is allocated for pilot.
% 16 point IFFT is taken and 4 samples of CP are added to generate samples corresponding to
% each time domain symbol.

%% Single carrier parameters
ovsf = 1;                         % Oversampling factor
fs = 168e3*ovsf;                  % Sampling rate
nfft = 16 * ovsf;                 % FFT size
cpLen = 4* ovsf;                  % CP length

scBin = 1;                        % Single subcarrier per subchannel
Nbins = 12;                       % Number of subchannels
Nsc = Nbins*scBin;                % Total number of subcarriers
lG = (nfft-Nsc)/2;                % number of left guard subcarriers
rG = lG;                          % number of right guard subcarriers
Nsym = 10000;                     % Number of time domain OFDM symbols
mod_data = '256QAM';              % 256QAM modulation for data
mod_pilot = 'BPSK';               % BPSK modulation for pilots
bitsPerSym = 8;                   % No of bits per symbol for 256QAM
bitsPerSym_pilot = 1;             % No of bits per symbol for BPSK
NAbins = 1;                       % Number of allocated subchannels
dataIdx = [lG+7];                 % Middle single subchannel occupied

%% Tx waveform generation
% The Tx signal is generated by subchannel mapping followed by IFFT and CP addition.
Nsym_pilot = floor(Nsym/9);      % No of pilot symbols
Nsym_data = Nsym - Nsym_pilot;   % No of data symbols

% Data bits
txBits = randi([0,1],1,NAbins*(scBin)*bitsPerSym*Nsym_data);
% Data symbols
dataSyms1 = qammod(txBits(:),256,'InputType','bit'); % 256QAM modulation
dataSyms(1,1:Nsym_data) = dataSyms1./sqrt(170); % normalising the symbols
scatterplot(dataSyms);title('Tx symbols');

% Pilot bits
txBits_pilot = randi([0,1],1,NAbins*(scBin)*bitsPerSym_pilot*Nsym_pilot);
% Pilot symbols
pilotSyms = (txBits_pilot - 1/2)*(-2); % BPSK modulation

txGrid = zeros(nfft,1);
txSig = [];
jj = 1; kk = 1;
for k1 = 1:Nsym
    symbol_index = mod(k1,9);
    if symbol_index == 5    % the middle symbol of 9 consecutive symbols is for pilot
        pilot_symbols = pilotSyms((jj-1)*NAbins*scBin+(1:NAbins*scBin));
        data_pilot_per_sym = pilot_symbols;
        jj = jj + 1;
    else
        data_symbols = dataSyms((kk-1)*NAbins*scBin+(1:NAbins*scBin));
        data_pilot_per_sym = data_symbols;
        kk = kk + 1;
    end
    % mapping to resource grid
    txGrid(dataIdx) = data_pilot_per_sym;
    % Tx waveform generation using IFFT
    txSym = ifft(fftshift(txGrid));
    % CP addition
    txSig = [txSig;txSym(nfft-cpLen+1:nfft); txSym];
    
end

%% FIR LPF for Tx Filtering
FilterOrder = 256;    % Number of taps-1
fpass = 7.05e3;       % Passband edge frequency in Hz
fstop = 10.09e3;      % Stopband edge frequency in Hz
Apass1 = 0.01;        % Passband ripple in dB
Astop1 = 90;          % Stopband attenuation in dB
firlp = dsp.LowpassFilter('SampleRate',fs, ...
    'FilterType','FIR', ...
    'DesignForMinimumOrder',false, ...
    'FilterOrder',FilterOrder, ...
    'PassbandFrequency',fpass, ...
    'StopbandFrequency',fstop, ...
    'PassbandRipple',Apass1, ...
    'StopbandAttenuation',Astop1);
C = coeffs(firlp);
firTaps = C.Numerator;
hFilter = firTaps;

%% FIR filtering of Tx signal
txSigFilt  = conv(txSig,hFilter);
txSigFilt  = txSigFilt(127:end); % delay correction introduced by 257 taps filter

figure;
SpectrumAnalyzer1(txSig, fs, 0, fs/1024, 'r');
hold on
SpectrumAnalyzer1(txSigFilt, fs, 0, fs/1024, 'g');grid on;
legend('Without filtering','With filtering','FontSize',11);
title('Tx Signal');

%% Plotting spectrum with FCC Mask
% Steps for plotting FCC Spectral Emission Mask
% Step 1. The spectrum of the signal is first plotted with 200 Hz resolution bandwidth.
% Step 2. Mark the line that is 26 dB below the maximum PSD value at 0 Hz from the above.
% Step 3. Obtain the bandwidth corresponding to the 26 dB line which is the
%    emission bandwidth.
% Step 4. Plot the spectrum with three times the emission bandwidth obtained in step 3.
% Step 5. The peak value of this PSD plot gives the peak power in dB of the
%    spectrum emission mask.
% Step 6. Plot the spectrum again with resolution bandwidth of 100 Hz and plot the mask
%    on top of it.

figure;
% Step 1
SpectrumAnalyzer1(txSigFilt, fs, 0, 200, 'b');grid on;
title('Tx Spectrum with 200 Hz res. BW');
% Step 2
[freq,spectrum]=SpectrumAnalyzer1(txSigFilt, fs, 0, 200, 'b');grid on;
freq0Idx = find(freq == 0);
spectrum0Hz = spectrum(freq0Idx);
spectrum26dBlow = (spectrum0Hz -26)*ones(1,length(freq));
hold on; plot(freq,spectrum26dBlow,'k');
legend('Single carrier with filtering','-26 dB line','FontSize',11);
% Step 3 : Obtain from the above plot
% Step 4
figure;SpectrumAnalyzer1(txSigFilt, fs, 0, 7.935*2*3*1e3, 'b');grid on;
title('Tx Spectrum with res. BW which is 3 times emission BW');
% Step 5 : Obtain from the above plot
% Step 6
figure;SpectrumAnalyzer1(txSigFilt, fs, 0, 100, 'b');grid on;
hold on; Mask_D(-11.2564);
legend('Spectrum with FIR filtering','FCC Mask D','FontSize',11);
title('Tx Spectrum with FCC Mask');

%% Noise addition (optional)
% Adding AWGN
SNR = 50;

% Noise setup
SNR_lin = 10^(SNR/20);    % Linear SNR

% Normalize noise power to take account of sampling rate, which is
% a function of the IFFT size used in OFDM modulation, and the
% number of antennas
N0 = 1/(sqrt(2.0*nfft)*SNR_lin);

% Create additive white Gaussian noise
noise = N0*complex(randn(size(txSigFilt)), ...
    randn(size(txSigFilt)));

% Add AWGN to the received time domain waveform
% RxSig = txSigFilt + noise ; % to add noise
RxSig = txSigFilt;

%% Recovery of Tx symbols : CP removal, applying FFT per symbol
offset = 0;
fftStart = cpLen;
RxSig_grid = [];

for ii = 1:Nsym
    RxSig_symbol = RxSig(offset+fftStart+(1:nfft));
    RxSig_fft_full = fftshift(fft(RxSig_symbol));
    RxSig_fft = RxSig_fft_full(dataIdx);
    offset = offset + nfft + cpLen;
    RxSig_grid = [RxSig_grid RxSig_fft];
end

%% Data symbol extraction
jj = 1;
for k1 = 1:Nsym
    symbol_index = mod(k1,9);
    if symbol_index ~= 5    % the middle symbol of 9 consecutive symbols is for pilot
        RxSymbols_data(:,jj) = RxSig_grid(:,k1);
        jj = jj+1;
    end
end

RxSymbols = RxSymbols_data(:);
scatterplot(RxSymbols); title('Rx symbols');

%% Decoding symbols to bits
RxSymbols2 = RxSymbols.*(sqrt(170)); % To account for the normalisation factor
Rxbits = qamdemod(RxSymbols2,256,'OutputType','bit');

%% SNR, EVM and BER calculation
% Measured SNR calculation
err_vec = RxSymbols - dataSyms(:);
noise_pow = sum(abs(err_vec).^2);
Rx_pow = sum(abs(RxSymbols).^2);
sig_pow = Rx_pow - noise_pow;
sigPowdB = 10*log10(sig_pow);
snr_db = 10*log10(sig_pow/noise_pow)

% EVM Calculation
normVal = sum(abs(dataSyms(:).^2));
evmNormalised = sqrt(noise_pow/normVal);
evmdB = 20*log10(evmNormalised)

% BER calculation
num_err = length(find(txBits(:) ~= Rxbits));
ber = num_err/(length(Rxbits))
