%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 Sigma Delta Modulator with DEM Algorithm              %%
%%                          Examples 8.3 & 8.4                           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t0=clock;

%-------------------------------------------------------------------------%
%                          Global variables                               %
%-------------------------------------------------------------------------%

R=64;                       % OSR
in_lvl =-30;                % input in dB
Fs=20e6;                    % Oversampling frequency
bw=Fs/2/R;
Ts=1/Fs;
N=1024*4;                   % Samples number
Vref=1;
Tstop=Ts*(300);
nper=11;
Fin=(nper*Fs/N);             % Input signal frequency (Fin = nper*Fs/N)
Ampl=Vref*20^(in_lvl/20);
Ntransient=0;
finrad=Fin*2*pi;             % Input signal frequency in radians
offset=0;

% *********************************************************************** %
%                           System parameters                             %
% *********************************************************************** %

VfsL=-1;
VfsH=1;

flag_quant=0;       %if 0 intervals equal
                    %if 1 first and last intervals are half LSB
                    
flag_dac=0.5;       % if 0  dac_output beginning of the interval
                    % if 0.5 dac_output half of the interval
                    % if 1   dac_output end of the interval
                    
DEM_sel=4;          % if 0 no algorithm
                    % if 1 DWA Algorithm
                    % if 2 ILA Algorithm (Rotation Approach)
                    % if 3 ILA Algorithm (Addition Approach)
                    % if 4 BUTTERFLY Algorithm
 
%%%%%%%%%%%%%%%%%%%%%%%
%  Flag nochange
%%%%%%%%%%%%%%%%%%%%%%%

nochange=0;             % change if the flag is 0

%%%%%%%%%%%%%%%%%
if nochange==0 
    clear Rnom err_Res Rreal vth vthidt vthid mism
NintQ=8;           % Number of intervals in the quantizer, counts two half quantizers equivalent to 1
Nelements=8;       % number of unary elements in the DAC

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Thresholds of the ADC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Delta=(VfsH-VfsL)/(NintQ-flag_quant);
N_Res=NintQ+flag_quant;
Rnom=ones(1,N_Res);
if flag_quant==1
    Rnom(1)=0.5;
    Rnom(N_Res)=0.5;
end
err_Res=0.00*randn(1,N_Res);
Rreal=Rnom.*(1+err_Res);
for i=1:1:N_Res-1
    vth(i)=VfsL+(VfsH-VfsL)*sum(Rreal(1:i))/sum(Rreal);
end

for i=1:1:N_Res-1
    vthidt(i)=VfsL+(VfsH-VfsL)*sum(Rnom(1:i))/sum(Rnom);
end
vthid=[VfsL vthidt VfsH]
if flag_quant==1
     vthid=[VfsL-Delta/2 vthidt VfsH+Delta/2];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unity elements of the DAC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mism=0.004*randn(1,Nelements);         % mismatches in the unity elements
mism=mism-sum(mism)/Nelements;
LSB=(VfsH-VfsL)/(NintQ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

if DEM_sel==2
      if Nelements>16
        disp('      ********************************************');
        disp('      **                                        **');
        disp('      **   Nelements must be lower than 16!!!   **');
        disp('      **                                        **');
        disp('      ********************************************');
        return
    end
end
if DEM_sel==3
      if Nelements>16
        disp('      ********************************************');
        disp('      **                                        **');
        disp('      **   Nelements must be lower than 16!!!   **');
        disp('      **                                        **');
        disp('      ********************************************');
        return
    end
end
if DEM_sel==4
    ctrl_butter=mod(log2(Nelements),1);
    if ctrl_butter~=0
        disp('      ********************************************');
        disp('      **                                        **');
        disp('      **   Nelements must be a power of 2!!!    **');
        disp('      **                                        **');
        disp('      ********************************************');
        return
    end
end


% ************************************************************************
%                       Open Simulink diagram first
% ************************************************************************

options=simset('RelTol', 1e-3, 'MaxStep', 1/Fs);
sim('SD2_DEM', (N+Ntransient)/Fs, options);	


% ************************************************************************
%   Calculates SNR and PSD of the bit-stream and of the signal
% ************************************************************************
w=hann(N);
f=Fin/Fs;			% Normalized signal frequency
fB=round(N*(bw/Fs));		% Base-band frequency bins
yy1=zeros(1,N);
yy1=yout(2+Ntransient:1+N+Ntransient)';
yy1=yy1-mean(yy1);
ptot=zeros(1,N);
[snr,ptot]=calcSNRdef(yy1(1:N),f,fB,w',N);
Rbit=(snr-1.76)/6.02;	% Equivalent resolution in bits


% ************************************************************************
%   Calculates Error Spectrum
% ************************************************************************
w=hann(N);
f=Fin/Fs;                % Normalized signal frequency
fB=round(N*(bw/Fs));     % Base-band frequency bins
yy2=zeros(1,N);
yy2=yout1(2+Ntransient:1+N+Ntransient)';

ptot1=zeros(1,N);
[snr1,ptot1]=calcSNRdef(yy2(1:N),f,fB,w',N);

% ************************************************************************
%  Graphic Output
% ************************************************************************
ptot=ptot-max(ptot);
ptot1=ptot1-max(ptot1);

figure(1);
semilogx(linspace(0,Fs/2,N/2), ptot1(1:N/2), 'r');
grid on;
title('PSD  of the Mismatch Error (DWA with 8 elements)')
xlabel('Frequency [Hz]')
ylabel('PSD [dB]')
axis([0 Fs/2 -140 -60]);
figure(4);
semilogx(linspace(0,Fs/2,N/2), ptot(1:N/2), 'b');
grid on;
title('PSD of the Output')
xlabel('Frequency [Hz]')
ylabel('PSD [dB]')
axis([0 Fs/2 -140 0]);


s1=sprintf('   In_level=%1.3f,    SNR(dB)=%1.3f,    Rbit=%1.3f', in_lvl, snr,  Rbit);
s2=sprintf('');
max(yout)
%s2=sprintf('   Simulation time =%1.3f min',etime(clock,t0)/60);
disp(s1)
disp(s2)
