%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% Matlab file for processing results of DSM1 in lib CT 
%% CTDSM Library
%% Mixed Signal IC Design Group                         
%% Boise State University, 2011.                          
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all; close all; clc;

 fs = 50e6;
 Nfft = 2^11;
 BW = fs/(2);     % Signal bandwidth
 tone_bin = max(primes(floor((BW/(5*fs))*Nfft))); % Find a tone bin
 fin = (tone_bin/Nfft)*fs;
 sim_time = Nfft*(1/fs)



% Define the result directory. Make sure thepath is correct and typo free.
% Using config simulation view. Switch to schematic if you simulate

rdir = '/tmp/kzhu/simulation/pipeline_top_sim/spectre/config/psf'; %% For spectre

% Pick transient data-set

dataset_name = 'tran-tran'; % For Spectre
% dataset_name = 'tran'; % For Ultrasim

% Extract all the signals

% signals = cds_srr(rdir, dataset_name);


% Extract required signals only% Extract all the signals
DAC_OUT = cds_srr(rdir, dataset_name, 'DAC_OUT'); vout = DAC_OUT.V; 
% Vinm = cds_srr(rdir, dataset_name, 'vinm'); vinm = Vinm.V; 
% vin = vinp - vinm;

%Clk = cds_srr(rdir, dataset_name, 'CLK'); clk = Clk.V;

%vout_dm = cds_srr(rdir, dataset_name, 'vout_dm'); vout = vout_dm.V;


%% Process data
VDD = 1.2;
VCM = VDD/2;
fmax = fs;  % Plot only till fmax
% N=length(vout);
 
N = Nfft;
vout = vout(1:N);
dt = 1/fs;
t1 = 0:dt:(N-1)*dt;  %time array
frange = 1/dt; % estimate the frequency range of the data

%% Apply Hann Window
%w = hann(N);
w = ones(N,1); % rectangular window
w1 = norm(w,1);
Y = fft(vout)/(N/2); %Apply window to estimate FFT
N1 = floor(N*fmax/frange); %Plot range (up to fmax)
Yl = Y(1:N1); %lower half of fft
DC1=Yl(1);  DC2 = mean(vout);  %Scaling test: DC1 should equal DC2

%% Estimate PSD
Py = (abs(Yl).^2)/(max(abs((Yl).^2))); %Normalized Power spectrum
max_bin = round(N1/1);
[Pmax, Pindex] = max(Py(1:max_bin)); %The tone power and index
fbin = Pindex-1;
fres = frange/N;  %freq step size and resolution
f = 0:fres:(N1-1)*fres;  %freq array
f1 = f/10^6 %change unit in x-axis

%% Estimate SNR and Neff
nb = 1; 
signal_bins = fbin + [-(nb-1)/2:(nb-1)/2];
inband_bins = 0:N1/(2*1);
noise_bins = setdiff(inband_bins,signal_bins);

SNR = 10*log10(sum(Py(signal_bins+1))/sum(Py(noise_bins+1)))
Neff = (SNR-1.76)/6.02;

%% Plot Waveform
if(1)
    figure()
    %plot(DAC_OUT.time, DAC_OUT.V)  % plot waveform
    plot(vout)  % plot waveform
    xlabel ('time'); ylabel ('DAC_out');
    title ('Pipeline ADC Output Waveform'); 
    grid on;
end; % if(0)

%% Plot Spectrum
if(1)
    figure()
    plot(f1, 10*log10(Py))  %plot spectrum
    xlabel ('Frequency (MHz)'); ylabel ('Normalized PSD (dB)');
    title ('Pipeline ADC recovered analog Output Spectrum'); 
    grid on;
    xlim([0 fs/2E06]);
end; % if(0)

    figure()
    semilogx(f, 10*log10(Py), 'r')  %plot spectrum
    text_handle = text(1e6,-40, ...
             sprintf('SQNR = %4.1f dB',SNR),...
        'FontWeight','bold','EdgeColor','red','LineWidth',3,'BackgroundColor',[1 1 1],'Margin',10);
    xlabel ('Frequency[Hz]'); ylabel ('dB');
    grid on;
   xlim([5e5 fs]);
   
   
break;
