% Select one of the Release 15 NR-TMs for FR1 and FR2 among:
% "NR-FR1-TM1.1","NR-FR1-TM1.2","NR-FR1-TM2",
% "NR-FR1-TM2a","NR-FR1-TM3.1","NR-FR1-TM3.1a",
% "NR-FR1-TM3.2","NR-FR1-TM3.3","NR-FR2-TM1.1",
% "NR-FR2-TM2","NR-FR2-TM2a","NR-FR2-TM3.1","NR-FR2-TM3.1a"

% or
% Select one of the Release 15 FRCs for FR1 and FR2 among:
% "DL-FRC-FR1-QPSK","DL-FRC-FR1-64QAM",
% "DL-FRC-FR1-256QAM","DL-FRC-FR2-QPSK",
% "DL-FRC-FR2-16QAM","DL-FRC-FR2-64QAM"

rc = "NR-FR1-TM3.2"; % Reference channel (NR-TM or FRC)

% Select the NR waveform parameters
bw = "10MHz"; % Channel bandwidth
scs = "30kHz"; % Subcarrier spacing
dm = "FDD"; % Duplexing mode

targetRNTIs = [];

displayEVM = true;
plotEVM = true;

if displayEVM
    fprintf('Reference Channel = %s\n', rc);
end

evm3GPP = false;

phaseNoiseOn = true;
IQImbalanceON = true;
filterOn = true;
nonLinearityModelOn = true;

OSR = 5; % oversampling factor

% Create waveform generator object
tmwavegen = hNRReferenceWaveformGenerator(rc,bw,scs,dm);

% Waveform bandwidth
bandwidth = tmwavegen.Config.ChannelBandwidth*1e6;

if OSR > 1
    % The |Config| property in |tmwavegen| specifies the configuration of
    % the standard-defined reference waveform. It is a read-only property.
    % To customize the waveform, make the |Config| property writable.
    tmwavegen = makeConfigWritable(tmwavegen);

    % Increase the waveform sample rate by multiplying the nominal sample
    % rate with |OSR|
    nominalSampleRate = getNominalSampleRate(tmwavegen.Config);
    tmwavegen.Config.SampleRate = nominalSampleRate*OSR;
else
    filterOn = false;
end

% Generate the waveform and get the waveform sample rate
[txWaveform,tmwaveinfo,resourcesinfo] = generateWaveform(tmwavegen,tmwavegen.Config.NumSubframes);
sr = tmwaveinfo.Info.SamplingRate; % waveform sample rate

txWaveform = txWaveform/max(abs(txWaveform),[],'all');

txWaveform = repmat(txWaveform,2,1);

if phaseNoiseOn
    % Carrier frequency
    if tmwavegen.Config.FrequencyRange == "FR1" % carrier frequency for FR1
        fc = 4e9;
    else % carrier frequency for FR2
        fc = 30e9;
    end

    % Calculate the phase noise level
    foffsetLog = (4:0.2:log10(sr/2)-0.001); % model offset from 10e3Hz to
                                            % almost sr/2. To avoid
                                            % aliasing, the sample rate
                                            % must be greater than twice
                                            % the largest value specified
                                            % by FrequencyOffset
    foffset = 10.^foffsetLog;    % linear frequency offset
    PN_dBc_Hz = hPhaseNoisePoleZeroModel(foffset,fc,'C');
    figure; semilogx(foffset,PN_dBc_Hz);
    xlabel('Frequency offset (Hz)');
    ylabel('dBc/Hz');
    title('Phase noise model'); grid on

    % Apply phase noise to the waveform
    pnoise = comm.PhaseNoise('FrequencyOffset',foffset,'Level',PN_dBc_Hz,'SampleRate',sr);
    pnoise.RandomStream = "mt19937ar with seed";
    rxWaveform = pnoise(txWaveform);
    release(pnoise);
else
    rxWaveform = txWaveform; %#ok<UNRCH>
end

if IQImbalanceON
    amplitudeImbalance = 0.2;
    phaseImbalance = 0.5;
    rxWaveform = iqimbal(rxWaveform,amplitudeImbalance,phaseImbalance);
end

if filterOn
    % Create low-pass filter object
    LPF = dsp.LowpassFilter('SampleRate',sr, ...
                            'FilterType','IIR', ...
                            'PassbandFrequency',sr/2-(sr/2*0.6), ...
                            'StopbandFrequency',sr/2-(sr/2*0.5), ...
                            'PassbandRipple',0.7, ...
                            'StopbandAttenuation',60);

    % Plot the magnitude response of the low-pass filter
    [h,w] = freqz(LPF);
    figure
    plot(w/pi,mag2db(abs(h)));
    axis('tight');
    grid;
    title('Magnitude Response of the LPF')
    xlabel('Normalized Frequency (x \pi rad/sample)');
    ylabel('Magnitude (dB)');

    % Filter the waveform
    rxWaveform = LPF(rxWaveform);
    release(LPF);
end

if nonLinearityModelOn
    % Generate Rapp model object
    rapp = comm.MemorylessNonlinearity('Method','Rapp model');
    rapp.Smoothness = 1.55;
    rapp.OutputSaturationLevel = 1;

    % Plot nonlinear characteristic
    plotNonLinearCharacteristic(rapp);

    % Apply nonlinearity
    rxWaveform = rapp(rxWaveform);
    release(rapp);
end

scope = spectrumAnalyzer('SampleRate',sr,...
    'ChannelNames',{'Before impairments','After impairments'},...
    'Title', 'Waveform before and after impairments');
scope([txWaveform,rxWaveform]);
release(scope);

if dm == "FDD"
    nFrames = 1;
else % TDD
    nFrames = 2;
end

rxWaveform(1:nFrames*tmwaveinfo.Info.SamplesPerSubframe*10,:) = [];

cfg = struct();
cfg.Evm3GPP = evm3GPP;
cfg.TargetRNTIs = targetRNTIs;
cfg.PlotEVM = plotEVM;
cfg.DisplayEVM = displayEVM;
cfg.IQImbalance = IQImbalanceON;

% Compute and display EVM measurements
[evmInfo,eqSym,refSym] = hNRDownlinkEVM(tmwavegen.Config,rxWaveform,cfg);

function plotNonLinearCharacteristic(memoryLessNonlinearity)
    % Plot the nonlinear characteristic of the power amplifier (PA) impairment
    % represented by the input parameter memoryLessNonlinearity, which is a
    % comm.MemorylessNonlinearity Communications Toolbox(TM) System object.

    % Input samples
    x = complex((1/sqrt(2))*(-1+2*rand(1000,1)),(1/sqrt(2))*(-1+2*rand(1000,1)));

    % Nonlinearity
    yRapp = memoryLessNonlinearity(x);

    % Release object to feed it a different number of samples
    release(memoryLessNonlinearity);

    % Plot characteristic
    figure;
    plot(10*log10(abs(x).^2),10*log10(abs(x).^2));
    hold on;
    grid on
    plot(10*log10(abs(x).^2),10*log10(abs(yRapp).^2),'.');
    xlabel('Input Power (dBW)');
    ylabel('Output Power (dBW)');
    title('Nonlinearity Impairment')
    legend('Linear characteristic', 'Rapp nonlinearity','Location','Northwest');
end

function [sampleRate] = getNominalSampleRate(cfgObj)
    % Obtain nominal sample rate for the input parameter CFGOBJ. CFGOBJ is
    % an object of type 'nrDLCarrierConfig'
    sampleRate = nr5g.internal.wavegen.maxSampleRate(cfgObj);
end