function cal_snr(N_samples, Fs, F_signal, ABCDs, nlev, umax, bandwidth, bin_ideal)

t = (0:N_samples-1) / Fs;
Amplitude = 0.75 * umax;  % Input at -2.5 dBFS to prevent loop saturation
u = Amplitude * sin(2 * pi * F_signal * t); 

% Execute Modulator Simulation using the scaled ABCD loop matrix
[v, xn, xmax, y] = simulateDSM(u, ABCDs, nlev);

%% 5. Windowed FFT & Spectral SNR Metric Extraction
% Apply a high-rejection Hann window to isolate bins accurately
win = hann(N_samples)';
v_windowed = (v - mean(v)) .* win; % Strip DC offset

% Run FFT and convert to power spectral scale
V_fft = fft(v_windowed);
V_mag = abs(V_fft(1:N_samples/2)) / (N_samples / 4);
V_power = V_mag.^2;

% Identify signal vs noise bins inside the target bandwidth
inband_bins = 1:floor((bandwidth / Fs) * N_samples);
signal_bin  = bin_ideal + 1; % Account for MATLAB 1-indexing
signal_span = 3;             % Number of bins signal power leaks into

% Isolate Signal Power
signal_bins_range = (signal_bin - signal_span) : (signal_bin + signal_span);
P_signal = sum(V_power(signal_bins_range));

% Isolate Noise Power (All in-band bins EXCEPT the signal ones)
noise_bins_range = setdiff(inband_bins, signal_bins_range);
P_noise = sum(V_power(noise_bins_range));

% Calculate SNR in decibels
SNR_db = 10 * log10(P_signal / P_noise);
fprintf('Matlab Calculated Modulator SNR for %d Hz signal : %.2f dB\n', F_signal, SNR_db);

end

