% clc;
% clear;
% close all;

%% =====================================================
% TEST SIGNAL GENERATION: QAM + RRC (FINAL VERSION)
% EXACTLY 1e5 SAMPLES
% =====================================================

%% PARAMETERS
FS        = 2e9;     % Sampling frequency (Hz)
M         = 16;        % QAM order

sps       = 20;        % Samples per symbol
rolloff   = 0.25;      % RRC roll-off factor
span      = 10;        % Filter span (symbols)

f_offset  = 0e6;       % Frequency offset (Hz)

target_samples = 1e5;  % FINAL REQUIRED LENGTH

rng(42);               % Reproducibility

%% =====================================================
% COMPUTE REQUIRED SYMBOLS
%% =====================================================
filter_delay = span * sps;  

% Over-generate to ensure enough samples
Nsym = ceil((target_samples + filter_delay) / sps);

fprintf('Using Nsym = %d\n', Nsym);

%% =====================================================
% PULSE SHAPING FILTER (RRC)
%% =====================================================
rrcFilter = rcosdesign(rolloff, span, sps, 'sqrt');

%% =====================================================
% DESIRED SIGNAL
%% =====================================================
data_x = randi([0 M-1], Nsym, 1);
x_sym = qammod(data_x, M, 'UnitAveragePower', true);

x = upfirdn(x_sym, rrcFilter, sps, 1);

%% =====================================================
% INTERFERER SIGNAL
%% =====================================================
data_i = randi([0 M-1], Nsym, 1);
i_sym = qammod(data_i, M, 'UnitAveragePower', true);

i_sig = upfirdn(i_sym, rrcFilter, sps, 1);

%% =====================================================
% REMOVE FILTER DELAY
%% =====================================================
delay = span * sps / 2;

x = x(delay+1:end-delay);
i_sig = i_sig(delay+1:end-delay);

%% =====================================================
% MATCH LENGTHS
%% =====================================================
min_len = min(length(x), length(i_sig));
x = x(1:min_len);
i_sig = i_sig(1:min_len);

%% =====================================================
% TRIM TO EXACT 1e5 SAMPLES
%% =====================================================
x = x(1:target_samples);
i_sig = i_sig(1:target_samples);

fprintf('Final signal length = %d samples\n', length(x));

%% =====================================================
% NORMALIZATION
%% =====================================================
x = x / rms(x);
i_sig = i_sig / rms(i_sig);

%% =====================================================
% FREQUENCY SHIFT
%% =====================================================
N = length(x);
t = (0:N-1)' / FS;

i_shifted = i_sig .* exp(1j * 2*pi * f_offset * t);

%% =====================================================
% PAPR (Crest Factor)
%% =====================================================
papr_x = 10*log10(max(abs(x).^2) / mean(abs(x).^2));
papr_i = 10*log10(max(abs(i_shifted).^2) / mean(abs(i_shifted).^2));

fprintf('PAPR (Desired)    = %.2f dB\n', papr_x);
fprintf('PAPR (Interferer) = %.2f dB\n', papr_i);

%% =====================================================
% SAVE SIGNALS
%% =====================================================
save('desired_signal.mat', 'x','FS');
save('interferer_signal.mat', 'i_shifted','FS');

%% =====================================================
% POWER SPECTRAL DENSITY (WELCH)
%% =====================================================
nfft = 4096;
window = hamming(1024);
noverlap = 512;

[p1, f1] = pwelch(x, window, noverlap, nfft, FS, 'centered'); 
[p2, f2] = pwelch(i_shifted, window, noverlap, nfft, FS, 'centered'); 

%% Desired Signal PSD
figure;
plot(f1/1e6, 10*log10(p1), 'LineWidth', 1.2);
grid on;
xlabel('Frequency (MHz)');
ylabel('PSD (dB/Hz)');
title('PSD of Desired Signal');

%% Interferer PSD
figure;
plot(f2/1e6, 10*log10(p2), 'LineWidth', 1.2);
grid on;
xlabel('Frequency (MHz)');
ylabel('PSD (dB/Hz)');
title('PSD of Interferer');

%% =====================================================
% EXPORT IQ DATA
%% =====================================================
x_export = x / max(abs(x));
i_export = i_shifted / max(abs(i_shifted));
