image thumbnail

Sigma Delta Test Bench

by

 

A test bench to simulate and analyse Sigma Delta modulators

Run_Me.m
clear; %clear the matlab environment..


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Modulator synthesize settings (For toolbox)

    %Modulator specification
    order = 1; % Order of the Modulator
    OSR = 64; % Oversampling ratio of the modulator (Used for toolbox  functions and behavioural simulation in simulink)
    form = 'CIFB'; %Modulator architecture. %At present, simulink simulations support only this structure. Other structures can be incorporated once 
    nLev = 2; % Number of levels in the quantiser 
    OBG = NaN; % Out of band gain of NTF (max 2 for binary quantiser), Use NaN for default by toolbox(OBG = 1.5)
    opt = 0; %Optimisation of NTF zeroes (0 = no opt, 1 = opt)
    
%Behavioural simulation settings (For Simulink)

%Inputs to the simulator
    Sinusoid_DC = 0; % 1 for Sinusoid input 0 for DC input 
    Fin_bin =35;
    Amplitude = 0.5; % Input sine wave amplitude or DC amplitude (Range 0 - 1) [ 1/127, 1/131, rect 8K and 22050/4 Sine
    %1/500 Idle Tones
    Fo = 22050; %Bandwidth of baseband in Hz
    NFFT = 8192; % Number of FFT points used for simulink data analysis (No windowing implemented at present, more points, better results)
    %As NFFT increases simulink simulation time increases (number of
    %samples gathered is will be greater than or equal to NFFT) [Schreier
    %suggests using atleast 64*OSR as NFFT
    Quantisation_Interval = 2/(nLev-1); %Determines the quantisation step size. Construct input signal based on this. 
    %Set quantisation interval param to 1 for 1 bit quantiser.
    random_num =[rand rand rand rand rand]; %DAC randomiser matrix  CTBD
    Windowed_FFT = 1; % 0 = FFTnot windowed, 1  = FFT windowed. Enable plots and disable stems since window compensation is done only for plots
    %If Fin is exactly on a bin, windowing is disabled CTBD
    Integrator_Gain = 1000000; % DC gain of integrator block 
    Dither_Power =0.00001; %Dither Gain (not in dB)
    Amplitude_scaling = 1; %(1 == Scaling) (0 == No Scaling) Scale amplitude to MSA reported by toolbox (For stable operation)
    
%Enable/Disable plots  (ON/OFF as required)
    Plot_Signals = 'OFF'; %ON = plot signals OFF = do not plot (Switch OFF plots while simulating low frequencies)
    Plot_Spectrum = 'ON'; % ON/OFF plotting spectrum
    Stem_plot = 0; % Use either Stem(1) or Plot(0) function to show spectrum of inputs for DC system input
    Plots_from_toolbox = 'OFF'; % ON/OFF plots from Toolbox (Plots == DB plots, Stem == Magnitude plots)
    

%Simulation runtime uncomment as needed
    
%Simulink Model to be used
if order == 1 && strcmp(form,'GEN')
    mdl = 'First_Order_Model_gen_new'; %File name of the simulink model without extension. Below mentioned parameters should be properly programed inside each model block for successful simulation
elseif order == 1 && strcmp(form,'CIFB')
    mdl = 'First_Order_Model'; %File name of the simulink model without extension. Below mentioned parameters should be properly programed inside each model block for successful simulation
elseif order == 2 && strcmp(form,'GEN')
    mdl = 'Second_Order_Model_Gen_new';    
elseif order == 2 && strcmp(form,'CIFB')
    mdl = 'Second_Order_Model';
elseif order == 3
    mdl = 'Third_Order_Model';
end


%Quantiser selection
if nLev == 2
    bin_quant = 1;
    multi_quant = 0;
elseif nLev > 2
    bin_quant = 0;
    multi_quant = 1;
else
    uiwait(msgbox('Invalid quantiser levels','Sigma Delta Simulator'));
    return;
end

%Parameter Verification
if order < 0
    uiwait(msgbox('Invalid modulator order','Sigma Delta Simulator'));
    return;
end

if OSR~= ceil(OSR)
    uiwait(msgbox('OSR should be an integer','Sigma Delta Simulator'));
    return;
end

if ~(strcmp(form,'CIFB') ||strcmp(form,'CRFB') ||strcmp(form, 'CIFF') ||strcmp(form,'CRFF') ||strcmp(form,'GEN')) 
    uiwait(msgbox('Unknown modulator form','Sigma Delta Simulator'));
    return;
end

if (Sinusoid_DC ~= 1)&&(Sinusoid_DC ~= 0)
    uiwait(msgbox('Only sinusoids and DC can be simulated','Sigma Delta Simulator'));
    return;
end

if (Fin_bin ~= ceil(Fin_bin))
    uiwait(msgbox('Please enter an integer for frequency bin','Sigma Delta Simulator'));
    return;
end

if (Amplitude >1) || (Amplitude < -1)
    uiwait(msgbox('Amplitude out of range','Sigma Delta Simulator'));
    return;
end

if ~((Amplitude_scaling == 1) || (Amplitude_scaling == 0))   
    uiwait(msgbox('Enter proper values for amplitude scaling','Sigma Delta Simulator'));
    return;
end    







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Parameters generated from inputs 
    Fsnyquist = 2*Fo; %Nyqusist sampling rate
    Tsnyquist = 1/Fsnyquist; %Nyquist sampling time.
    Fs = Fsnyquist*OSR; %Oversampling rate
    Ts = 1/Fs; % Sampling Interval
    Fin = Fs*Fin_bin/NFFT;
    Tcont = 1/(Fs*1000); % For plotting continous time signals this interval will be used along with variable time_index_cont
    Trun = 1/Fin; % Trun = 1 cycle of input frequency
    Nrun_sinusoid = (NFFT*2)*(Fin/Fs)+1; % Number of sinusoid cycles should be able to give atleast  NFFT*2 samples and Fs/Fin indicates number of samples per input cycle
    Trun_simulink = (NFFT*2)*Ts; % Simulink runs for Nrun cycles of input frequency
    Trun_simulink_DC = Nrun_sinusoid*Trun; % Simulink run time for DC inputs to gather NFFT*2 samples.
    NFFT_DEC = NFFT/OSR; % Number of FFT points for Decimated outputs
    fspacing = Fs/NFFT;

    %Trun = 1; % Trun in seconds
    
% Window initialization for FFT
    Window = hann(NFFT); %Hann window
    no_valid_signal_bins = 3; %Number of valid signal bins 
    w1 = norm(Window,1);
    w2 = norm(Window,2);
    
%String cocatenation for titile    
    title_str = strcat(strcat(' OSR: ',num2str(OSR)),strcat(strcat(', Order:',num2str(order)), strcat(' ,Quantiser Levels:',num2str(nLev)) )); %Concatenating strings for title      NBW = (w2/w1)^2; % Noise Bandwidth
    

%Code for Nth order Sig Delta Modulator
if ~strcmp(form,'GEN')
%Step 1, synthesize the NTF


%Step 1: Synthesize NTF
    %See order and OSR defentions at variables section on Top
    H = synthesizeNTF(order,OSR,opt,OBG);  % Synthesize the NTF for the modulator with order, OSR define
    
    
%Step 2: Realize NTF    
    %Realise the NTF with a particular
    %Modulator structure. Examples are CIFB, CIFF, CRFB, CRFF. 
    %We will be using a CIFB structure to establish the whole design path and
    %change if required later.
    [a,g,b,c] = realizeNTF(H,form); %This function converts the NTF synthesized earlier into coeffiecients for circuit implementation.
    
%Step 3: Generate loop filter description in ABCD format    
    %realizeNTF returns coefficients for an unscaled modulator whose internal
    %stages occupy an unspecified range. So dynamic range scaling must be
    %performed (Page 268 for a diagram)
    b(2:end) = 0; % Set all feed in terms to 0 for a maximally flat STF (??)
    ABCD = stuffABCD(a,g,b,c,form);  % ABCD matrix (generic loop filter description) (unscaled)
    
%Step 4: Dynamic range scaling of ABCD matrix    
    %Use scaleABCD to determine scaling factors for each state of
    %modulator. Modulator is simulated with inputs of different amplitudes
    %to determine maximum stable input amplitude (umax). Then dynamic range
    %scaling is peformed to maximum value of each state equals specified
    %limit (xlim)
    xLim = 0.9;
    f0 = 0; %Normalised frequency of the test sinusoid
    %ABCDs is the scaled state space matrix description
    %umax is the maximum possible amplitude of the signal that can be
    %applied to the modulator.
    
    %Use either of one
    %[ABCDs umax] = scaleABCD(ABCD,nLev,f0,xLim);
    [ABCDs umax] = scaleABCD(ABCD);
    %If Amplitude scaling is enabled, user defined amplitude is scaled to
    %MSA levele (umax) reported by toolbox
    if ((Amplitude_scaling == 1) && (Amplitude > umax))
        Amplitude = umax*Amplitude;
    end
    
%Step 5: Convert scaled ABCD matrix back to coeffcients for circuit implementation     
    %Scaled ABCD matrix is translated back into co-efficients by mapABCD
    %for the modulator structure defined.
    [a g b c] = mapABCD(ABCDs,form);
    %This data goes into simulink models
    
% Additional toolbox functions used below for plotting and analysis    
    
    
    if strcmp(Plots_from_toolbox,'ON')
        fig0 = figure('NumberTitle','off','Name',strcat('Plots from toolbox  ',title_str));
        scrn = get(0,'ScreenSize');
        set(fig0,'Position',[scrn(3)/8 scrn(4)/16 0.75*scrn(3) 0.80*scrn(4)]);    
        clf;

        subplot(2,2,1);
        plotPZ(H);  %Plots the poles and Zeroes of NTF synthesized above
        title('Poles and Zeroes of NTF');
    end
    
    RMS_Gain_from_toolbox_in_band_of_interest = dbv(rmsGain(H,0,0.5/OSR)); % Calculates rms gain of NTF in the band of interest (0 to f0, f0 = 0.5/OSR)



    %Simulate the Delta Sigma Modulator with the NTF sythesized
   
    
    %Tone bin calculation. 
    % Bin  1        Bin 4096        Bin 8192        --> FFT bins (Example
    % showing with NFFT = 8192)
    % 0 Hz          11025 Hz        2822400 Hz      --> 0 to Fs (2*Fo*OSR)
    %0              0.5             1               --> Normalised frequency 
    % Bin - Bin frequency difference = Fs/8192
    bin_bin_freq_diff = Fs/NFFT; 
    
    %Bin number of input signal = Input frequency / bin-bin frequency
    %difference
    
    tone_bin = Fin/bin_bin_freq_diff;
    
    t = [0:NFFT-1]; %Timescale
    
    u = (umax*Amplitude)*(nLev-1)*sin(2*pi*(tone_bin/NFFT)*t); %Input signal generated  (Freq = tone_bin/NFFT)
    %Note that input signal amplitude is scaled to umax to prevent
    %modualtor becoming unstable
    
    %Frequency band of interest is 0 to F0, where F0 = Fs/2*OSR =
    %1/2*OSR=0.5/OSR = 0.0078. Selected input is withing this band.
    v=simulateDSM(u,H,nLev); % Simulates the DSM
    %To plot simulated results
    n = 1:1:Fs/Fin; %Index for samples to be displayed (Fs/Fin indicates samples required to show one period of input)
   
    if (length(n)>length(t)) %To enusre that plots below are not out of index. One complete cycle may not be displayed if n> time index
        n=1:length(t);
    end
    
    
    if strcmp(Plots_from_toolbox,'ON')
        subplot(2,2,2);
        stairs(t(n),v(n),'g'); %Input plotted against time
        grid on;
        hold on;
        stairs(t(n),u(n),'b');% Output plotted against time
        title('Input and Output');
        axis tight;
        xlabel('Time');
        ylabel('Amplitude');
        leg1 = legend('Modulator output','Input signal');
        set(leg1,'Location', 'SouthWest');
    end
    
    %Note Windowed FFT
    spec = fft(v.*hann(NFFT))/(NFFT*(nLev-1)/4); % FFT of the windowed output (Normalised )
    snr_toolbox = calculateSNR(spec(1:ceil(NFFT/(2*OSR))+1),tone_bin); % snr calculation from FFT for the particular frequency tone_bin
    %Note that tone_bin is rounded , if input frequency is not in an
    %integer bin, snr may not be correct.
    
    NBW = 1.5/NFFT; % Noise Bandwidth for 
    f = linspace(0,0.5,NFFT/2+1); % Frequency range to plot from 0 to 0.5Fs and NFFT/2 = 8192/2 = 4096 plus one more point = 4097 points
    Sqq = 4* (evalTF(H,exp(2i*pi*f))/(nLev-1)).^2/3; % TBD
    if strcmp(Plots_from_toolbox,'ON')
        subplot(2,2,3);
        plot(f, dbv(spec(1:NFFT/2+1)),'b'); %Plots FFT wrt Freq (This is from simulation of the designed NTF)
        grid on;
        hold on;
        plot(f, dbp(Sqq*NBW),'m', 'Linewidth',1); % Plots the expected PSD
        title(' Plots of  expected PSD and FFT of output');
        xlabel('Frequency');
        ylabel('Amplitude/Power (dB)');
        leg2 = legend('FFT of output','Expected PSD');
        set(leg2,'Location', 'SouthEast');
    end


    %Usage of simulate SNR function to give sinusoids of different amplitude to
    %and same frequency to the modulator to see the resules
    %TBD formulate this section to simulate from 1Hz to 20KHz and different
    %amplitudes and save the plots.x 
    
    %amp = [-130:5:-20 -17:2:-1]; %Amplitude range
    %snr = simulateSNR(H,OSR,amp,[],nLev); % SNR simulation for amplitudes listed above
    %figure(5);
    %clf;
    %plot(amp,snr,'-b',amp,snr,'db');  %SNR plot wrt amplitudes
    %title('SNR with different input amplitudes');
    %[pk_snr pk_amp] = peakSNR(snr,amp); 



    %To plot STF and NTF
    [Ha Ga] = calculateTF(ABCD); % Calculates STF from ABCD matrix. Ha = NTF eq, Ga = STF eq
    f = linspace(0,0.5,100);
    z = exp(2i*pi*f);
    magHa = dbv(evalTF(Ha,z)); 
    magGa = dbv(evalTF(Ga,z));
    
    if strcmp(Plots_from_toolbox,'ON')
        subplot(2,2,4);
        plot(f,magHa,'b',f,magGa,'m','Linewidth',1);
        grid on;
        title('STF and NTF')
        leg3 = legend('NTF','STF');
        set(leg3,'Location', 'SouthEast');
    end
   
end

uiwait(msgbox('Continue with behavioural simulation','Sigma Delta Simulator'));



%Behavioural Simulation using simulink
   
 
%Outputs from the system
    sampled_input = []; %Samples from the input to the modulator
    sampled_modulator_output = []; %Samples from Sigma Delta Modulator
    sampled_filter_output = []; %Samples from digital LPF
    sampled_decimated_output = []; %Decimated out put

if Sinusoid_DC == 1       

%Loading and starting simulink simulation 
    open_system(mdl); %Open simulink block
    set_param(mdl,'StartTime','0'); %Setting simulink simulation start time to 0 Seconds
    set_param(mdl,'StopTime',num2str(Trun_simulink)); %Setting simulink stop time  (Converted run time Trun to string as required).
    open_system(strcat(mdl,'/Scope2'));
    sim(mdl); % Run simulink simulation
    modulator_output_samples = sampled_filter_output; % Copy samples from simulink simulation results
    time_index = linspace(0,Trun_simulink,length(sampled_modulator_output)); 
    time_index_OSR = linspace(0,Trun_simulink,length(sampled_filter_output));
    time_index_input = linspace(0,Trun_simulink,length(sampled_input)); %Time index for plotting continous time signal
    time_index_decimated = linspace(0,Trun_simulink,length(sampled_decimated_output)); %Time index for plotting Decimated output.
    
elseif Sinusoid_DC == 0
    
%Loading and starting simulink simulation 
    open_system(mdl); %Open simulink block
    set_param(mdl,'StartTime','0'); %Setting simulink simulation start time to 0 Seconds
    set_param(mdl,'StopTime',num2str(Trun_simulink_DC)); %Setting simulink stop time  (Converted run time Trun to string as required).
    open_system(strcat(mdl,'/Scope2'));
    sim(mdl); % Run simulink simulation
    modulator_output_samples = sampled_filter_output; % Copy samples from simulink simulation results
    time_index = linspace(0,Trun_simulink_DC,length(sampled_modulator_output)); 
    time_index_OSR = linspace(0,Trun_simulink_DC,length(sampled_filter_output));
    time_index_input = linspace(0,Trun_simulink_DC,length(sampled_input)); %Time index for plotting continous time signal
    time_index_decimated = linspace(0,Trun_simulink_DC,length(sampled_decimated_output)); %Time index for plotting Decimated output.    
    
end
  
    
%Plotting input, intermediate signals and outputs    
if strcmp(Plot_Signals,'ON')
    fig1 = figure('NumberTitle','off','Name',strcat('Input,Modulator o/p, filter o/p and decimator o/p {',strcat(title_str, ' }')));
    scrn = get(0,'ScreenSize');
    set(fig1,'Position',[scrn(3)/8 scrn(4)/16 0.75*scrn(3) 0.80*scrn(4)]);
    clf;
    
    subplot(4,1,1);
    plot(time_index_input,sampled_input,'g','LineWidth',1); %Plotting system input frequency (Continous time domain) 
    grid on;
    axis tight;
    title('System Input (Sampled at F0*2*OSR)');
    %xlabel('Time');
    ylabel('Amplitude');  
    
    subplot(4,1,2);
    plot(time_index,sampled_modulator_output); %Plotting modulator output
    grid on;
    axis tight;
    title('Modulator Output');
    %xlabel('Time');
    ylabel('Amplitude');
    
    subplot(4,1,3);
    stairs(time_index_OSR,sampled_filter_output,'r','LineWidth',2); % Plotting Digital LPF output 
    grid on;
    axis tight;
    title('Digital Filter Output');
    %xlabel('Time');
    ylabel('Amplitude'); 
    
    subplot(4,1,4);
    stairs(time_index_decimated,sampled_decimated_output,'k'); % Plotting decimated output
    grid on;
    axis tight;
    title('Decimated Output');
    xlabel('Time');
    ylabel('Amplitude');
    %leg = legend('Modulator Output','Digital filter output','Modulator Input');
    %set(leg,'Location','SouthWest');
end
    
 %Spectrum analysis
 sampled_input_for_fft = sampled_input(NFFT/4:(NFFT/4)+(NFFT-1));
 Windowed_input = reshape(sampled_input_for_fft,1,length(sampled_input_for_fft)).*Window; %Windowing the input also reshaping input matrix to match window matrix

 
 sampled_output_for_fft = sampled_modulator_output(NFFT/4:(NFFT/4)+(NFFT-1));
 Windowed_output = Window.*reshape(sampled_output_for_fft,1,length(sampled_output_for_fft)); %Windowing the output also reshaping input matrix to match window matrix
 
 
 decimated_output_for_fft = sampled_decimated_output(NFFT_DEC/2:end);
 Window_decimated = hann(length(decimated_output_for_fft)); %Hann window
 Windowed_decimated_output = Window_decimated.*reshape(decimated_output_for_fft,1,length(decimated_output_for_fft)); %Windowing the decimated output also reshaping input matrix to match window matrix
 
 %FFT calculation

 input_spectrum = fft((Windowed_input)/(w1/2),NFFT); % FFT calculation
 
 modulator_spectrum = fft(Windowed_output/(w1/2),NFFT); % FFT calculation
 
 decimated_spectrum = fft(decimated_output_for_fft,NFFT_DEC); % FFT calculation without windowing
 
 
 
 for i = 1:NFFT % This block ensures that FFT results if 0, will be machine epsilon,eps  (later useful when converts to dB, else dB will give -Infinity
     if input_spectrum(i) == 0
         input_spectrum(i) = rand*eps; %Multiplied by rand only for visual effect of noise floor at ~ -300dB
     end
     if modulator_spectrum(i) == 0
         modulator_spectrum(i) = rand*eps; %Multiplied by rand only for visual effect of noise floor at ~ -300dB
     end
 end
 
 for i = 1:length(decimated_spectrum) % This block ensures that FFT results if 0, will be eps (later useful when converts to dB, else dB will give -Infinity
     if decimated_spectrum(i) == 0
         decimated_spectrum(i) = rand*eps; %Multiplied by rand only for visual effect of noise floor at ~ -300dB
     end
 end
 
  
  %Specturm of Inputs
  if strcmp(Plot_Spectrum,'ON')      

    %Frequency_axis_calculation 
    fspacing = Fs/NFFT; 
    fnyquist = Fs/2;
    fstart = -fnyquist;
    fend = fnyquist - fspacing;
    faxis = fstart:fspacing:fend;
    faxis_half = faxis((NFFT/2)+1:end);
    faxis_OSR = faxis_half(1:((NFFT/(2*OSR))+1)); % For frequencies upto F0
    
    faxis_log = [1:NFFT/2]/NFFT;
        
       
    fig2 = figure('NumberTitle','off','Name',strcat('Frequency Spectrum: ',title_str));
    scrn = get(0,'ScreenSize');
    set(fig2,'Position',[scrn(3)/16 scrn(4)/16 0.9*scrn(3) 0.80*scrn(4)]);    
    
   %Spectrum of modulator input upto F0 (Band of interest) 
    subplot(2,2,1);
    if Stem_plot == 0        
        plot(faxis_OSR,dbv(input_spectrum(1:NFFT/(2*OSR)+1)));
        ylabel('Amplitude (dB)');  
    elseif Stem_plot == 1
        stem(faxis_OSR,(input_spectrum(1:NFFT/(2*OSR)+1)));
        ylabel('Amplitude');  
    end
    grid on;
    axis tight;
    title('Single sided amplitude spectrum of Input upto F0');
    xlabel('Frequency (Hz)');

    
   %Spectrum of modulator input upto Fs/2 (F0*2*OSR / 2) 
    subplot(2,2,2);
    semilogx(faxis_log,dbv(input_spectrum(2:NFFT/2+1)));
    hold on;
    temp = log10([(NFFT/(2*OSR))+2]/NFFT);
    plot([temp temp],[min(dbv(input_spectrum(2:NFFT/2+1))) 0],'r','Linewidth',1.5);
    grid on;
    axis tight;
    ylabel('Amplitude (dB)');
    title('Single sided amplitude spectrum of Input upto Fs/2');
    xlabel('Normalised Frequency');
   
    
%Spectrum of outputs
   
 %Spectrum of modulator output upto F0
    subplot(2,2,3);
    if Stem_plot == 0
        plot(faxis_OSR,dbv(modulator_spectrum(1:NFFT/(2*OSR)+1))); 
        ylim([-120 0]);%Axis scaling
        ylabel('Amplitude (dB)');        
    elseif Stem_plot == 1
        stem(faxis_OSR,(modulator_spectrum(1:NFFT/(2*OSR)+1))); % Was working here change this for input
        ylabel('Amplitude');        
    end
    grid on;
    axis tight;
    title('Single sided amplitude spectrum of Modulator output upto F0');
    xlabel('Frequency (Hz)');
      
  %Spectrum of modulator output upto Fs/2 
    subplot(2,2,4);
    semilogx(faxis_log,dbv(modulator_spectrum(2:NFFT/2+1)));
    %hold on
    %[f p] = logsmooth(modulator_spectrum,Fin_bin,2,no_valid_signal_bins);
    %plot(faxis_half,p,'m','Linewidth',1.5);
    grid on;
    %ylim([-120 0]);%Axis scaling
    axis tight;
    ylabel('Amplitude (dB)');
    title('Single sided amplitude spectrum of Modulator output upto Fs/2');
    xlabel('Normalised Frequency');
        
  %Spectrum of Decimated outputs upto F0 

    decimated_spectrum_shift = fftshift(abs(decimated_spectrum))/NFFT_DEC;
    
    %Frequency axis calculation
    fspacing = Fsnyquist/NFFT_DEC;
    fnyquist = Fsnyquist/2;
    fstart = -fnyquist;
    fend = fnyquist - fspacing;
    faxis = fstart:fspacing:fend;
    faxis_half = faxis((NFFT_DEC/2)+1:end);
    
    decimated_spectrum_half = 2*decimated_spectrum_shift((NFFT_DEC/2)+1:end);%Spectrum amplitude multiplied by 2 since taking only single sided spectrum
    decimated_spectrum_half(1) = decimated_spectrum_half(1)/2; % For DC first FFT bin should not be multiplied by 2    
    
    
    %Spectrum plot
    figure('NumberTitle','off','Name',strcat('Frequency Spectrum (Decimated outputs): ',title_str));
    if Stem_plot == 0
        plot(faxis_half,mag2db(decimated_spectrum_half));
        %ylim([-120 0]);%Axis scaling
        ylabel('Amplitude (dB)');
    elseif Stem_plot == 1
        stem(faxis_half,(decimated_spectrum_half)); 
        ylabel('Amplitude');
    end
    grid on;
    axis tight;
    title('Single sided amplitude spectrum of Decimated output');
    xlabel('Frequency (Hz)');
  
end
 

%SNDR, ENOB and SFDR from decimated spectrum.
if (Sinusoid_DC == 1)
    %New SNDR
    signal_bins  = Fin_bin +[-(no_valid_signal_bins-1)/2:(no_valid_signal_bins-1)/2]; %Signal bins
    inband_bins  = 0:NFFT/(2*OSR); % All inband bins
    noise_bins   = setdiff(inband_bins,signal_bins);
    signal_power = sum(abs(modulator_spectrum(signal_bins+1)).^2);
    noise_power  = sum(abs(modulator_spectrum(noise_bins+1)).^2);
    
    SNDR     = dbp(signal_power/noise_power);
    SNDR_dBFS = dbp(1/noise_power);
    SFDR_dBFS     = -1*max(dbp(abs(modulator_spectrum(noise_bins+1)).^2));
    ENOB = (SNDR + 1.76)/6.02;
end


%Printing out results in matlab command window
if ~strcmp(form,'GEN')
    Max_Stable_Amplitude = umax
    RMS_Gain_from_toolbox_in_band_of_interest
end

if (Sinusoid_DC == 1)
    SNDR_from_simulink = SNDR
    SNDR_dBFS_simulink = SNDR_dBFS
    SFDR_dBFS_from_simulink = SFDR_dBFS
    ENOB_from_simulink_results = ENOB
    snr_toolbox
end
  
  
%NTF plots from toolbox and simulink merged

if ~(strcmp(form,'GEN') ||(Sinusoid_DC == 0))
    fig5 = figure('NumberTitle','off','Name','Comparison of outputs from toolbox and simulink');
    scrn = get(0,'ScreenSize');
    set(fig2,'Position',[scrn(3)/16 scrn(4)/16 0.9*scrn(3) 0.80*scrn(4)]);    
    
    %NTF comparison
    %f = logspace(0,0.5,NFFT/2+1);
    faxis_log = [1:NFFT/2]/NFFT;
    semilogx(faxis_log, dbv(spec(1:NFFT/2)),'b'); % NTF plot from toolbox
    hold on;
    semilogx(faxis_log,dbv(modulator_spectrum(2:NFFT/2+1)), 'm');
    grid on;
    ylim([-120 0]);%Axis scaling
    title('NTF comparison');
    xlabel('Frequency bins Normalised to Fs');
    ylabel('Amplitude (dB))');
    leg1 = legend('NTF (Toolbox)','NTF (Simulink)');
    set(leg1,'Location', 'SouthEast');
end    
    
    

Contact us