clc; clear all; close all;

%% System Parameters
fc = 1442.5e6;                             % Center frequency
numRBGs = 4;                               % Number of RBGs
rbgSize = 8;                               % RBs per RBG
modulation = 'QPSK';                      % Modulation scheme
numSubframes = 1;                          % Number of subframes
filenamePrefix = 'mimo_baseband';

a = -0.276;
b = 0.317; 
x = linspace(a, b, 10);
snum = 10; 
gain_dB_list = [-20, -10, 0, 10, 20] + x(snum);  % Gain levels in dB

%% 5G NR Downlink Configuration
cfgDL = nrDLCarrierConfig;
cfgDL.ChannelBandwidth = 50;              % MHz
cfgDL.FrequencyRange = 'FR1';
cfgDL.NumSubframes = numSubframes;
cfgDL.CarrierFrequency = fc;

cfgDL.SCSCarriers{1}.SubcarrierSpacing = 15;
cfgDL.SCSCarriers{1}.NSizeGrid = 150;
cfgDL.BandwidthParts{1}.SubcarrierSpacing = 15;
cfgDL.BandwidthParts{1}.NSizeBWP = 150;

% PRB allocation
rbIndices = 0:149;

%% Loop over gains
nfft = 4096;
for g = 1:length(gain_dB_list)
    gain_dB = gain_dB_list(g);
    gain_linear = 10^(gain_dB/20);

    cfgDL.PDSCH = {}; 
    pdsch = nrWavegenPDSCHConfig;
    pdsch.Modulation = modulation;
    pdsch.NumLayers = 1;
    pdsch.PRBSet = rbIndices; 
    pdsch.SymbolAllocation = [0, 14];
    cfgDL.PDSCH{1} = pdsch;
    
    % Generate waveform
    [waveform, info] = nrWaveformGenerator(cfgDL);
    waveform_scaled = gain_linear * waveform;
    numSamples = size(waveform_scaled, 1);
    fs = numSamples / (numSubframes * 1e-3);

    % Plot PSD
    [psd, f] = pwelch(waveform_scaled, hamming(nfft), [], nfft, fs, 'centered');
    f_abs = f + fc;
    
    % figure;
    % plot(f_abs/1e6, 10*log10(psd));
    % xlabel('Frequency (MHz)'); 
    % ylabel('PSD (dB/Hz)');
    % title(sprintf('PSD | Gain = %d dB', gain_dB)); 
    % grid on;
    % xline(1403, '--b', 'L-band Start');
    % xline(1430, '--b', 'L-band End');


    % Compute L-band Leakage
    bandMask = (f_abs >= 1400e6) & (f_abs <= 1427e6);
    f_band = f_abs(bandMask);
    psd_band = psd(bandMask);
    df = mean(diff(f_band));
    P_watt = sum(psd_band) * df;
    P_dBm = 10 * log10(P_watt / 1e-3);
    wl = physconst('Lightspeed')/fc; 


    Pr_dbm = P_dBm + 5 + 12.84 + 20*log10(wl/(4*pi));
    Prfi = 10 ^ ((Pr_dbm - 30)/10);

    Tb = (Prfi) / (physconst('Boltzmann') * 27 * 1e6); 
    fprintf('Gain = %.2f dB | Leaked Power into 1400–1427 MHz Band: %.2f dBm | Recieved Power: %.2f dBm\n', gain_dB, P_dBm, Pr_dbm);

    % Plot spectrogram
    window = 512; noverlap = 256; nfft_spec = 1024;
    [S,F,T,P] = spectrogram(waveform_scaled, window, noverlap, nfft_spec, fs, 'centered');
    
    % Plot
    figure;
    f_spec = F + fc; 
    freqMask = (f_spec >= 1400e6) & (f_spec <= 1427e6); 
    F_plot = f_spec(freqMask); 
    P_plot = P(freqMask); 
    imagesc(T*1e6, (F_plot)/1e6, 10*log10(abs(P_plot)));
    % imagesc(T*1e6, (F+fc)/1e6, 10*log10(abs(P)));
    
    % freqMask = (F >= -22e6) & (F <= 22e6); 
    % F_plot = F(freqMask); 
    % P_plot = P(freqMask); 
    % imagesc(T*1e6, (F_plot)/1e6, 10*log10(abs(P_plot)));

    axis xy;
    colormap jet;
    colorbar;
    clim([-170 -85]);
    title(sprintf('Spectrogram | Gain = %d dB', round(gain_dB)));
    xlabel('Time (\mus)'); 
    ylabel('Frequency (MHz)');
    % yline(1400, '-w', 'L-band Start');
    % yline(1427, '-w', 'L-band End');
    filename = sprintf('out_fc2_4RB_Gain%d_sn%d.png',round(gain_dB),snum);
    saveas(gcf, filename);

end
