image thumbnail

Data Weighted Averaging for Simulink

by

 

23 Feb 2009 (Updated )

A group of Delta-Sigma SIMULINK- models with DWA mismatch shaping.

run_dwa.m

% suggestions:
% [BPmode,DL,tf_mode] = [0 1 1] for lowpass
% [BPmode,DL,tf_mode] = [1 4 1] for bandpass (4-path DWA)
% [BPmode,DL,tf_mode] = [1 2 0] for bandpass (back-and-forth pointer)

Bits = 3;
BPmode = 0;     %stimulus and delta-sigma converter type
                % 0: lowpass,  1: bandpass
DACmode = 1;    % 1 FOR DELTA-SIGMA DAC, 
                % 0 FOR DELTA-SIGMA ADC
                
                
%---- DWA SETUP ----
DL = 1;         % delay in DWA
tf_mode=1;      % Determines the type of mismatch noise transfer function
                % 0: 1+z^(-DL),  1: 1-z^(-DL)
IDWAmode = 1;   % 0: normal DWA
                % 1: Incremental DWA (IDWA), i.e. ONE EXTRA UNIT ELEMENT

                
%DNL error related
Estd=0.01;      % standard deviation of DNL errors
shape = 'sine'; k = 3;  %the DNL error shape
% shape:    k
% 'poly'    polynomial order
% 'rand'    random error, k is negligible, e.g. 0
% 'sine'    sinusoid DNL error shape with k periods

N=pow2(Bits)-1;  %number of elems
if IDWAmode
    N=N+1;
end

random_dir = 1;  % pointer direction can change as a random pattern, as 
                 % long as the pattern length is D. See how the controlling
                 % sequence seq is created. If 1: random seq, else: ones.
%---- DWA SETUP END ----

if tf_mode==1
    Bf = [1,zeros(1,DL-1),-1];
else
    Bf = [1,zeros(1,DL-1),1];
end

%---- Simulation time
NT=pow2(17); %this many samples in simulation
fs = 1;ts= 1/fs; % sampling freq & interval
t=0:(NT-1)*ts;
Tfinal = t(end);
%----Simulation time END

%---- Stimulus
amp = 0.01;
fmult = 1/5;    % Input sinusoid's place in the signal band: 
                 % left side to right side 0->1. 0.5 is in the middle.
fb=fs*0.5/32;      % signal band edge in lowpass mode

fsig = fmult*fb; % signal freq
J=fsig*NT/fs;

if abs(J-round(J))>0
    J=round(J);
    if rem(J,2)==0
        J=J+1;
    end
    J = round(J);
end

if BPmode
    J=J+NT/4;
end

fsig = J/NT*fs;

data = amp*cos(2*pi*fsig*t);
%---- Stimulus END

%---- DELTA SIGMA ----
%Oversampling ratio is set to 32, order is 2 for lowpass, 4 for bandpass. 
%Also, loopfilter optimized for Bits = 3.
%Get DStoolbox to customize the DS loopfilter

if BPmode
        A = [-1.056673642799789  -2.056673642799789 0 0;1 1 0 0;1 1 -0.943326357200210  -1.943326357200210;0 0 1 1];
        B = [1 -1;0 0;0 0;0 0];
        C = [0  -0.612223827873632  -1.897128681596029  -1.931825636131040];
        D = [0 0];
else
        A = [0.996788098211802  -0.003211901788198;1   1];
        B = [1    -1;0     0];
        C = [2.446130189023710   1.771892689428995];
        D = [0 0];
end
%---- DELTA SIGMA END----

fsig = J/NT;

WLscale=0.5*(pow2(Bits)-1);  %
WLoffset = sum(0:(pow2(Bits)-1))/pow2(Bits);

K=0:(N-1);
simin=[t.' data'];

%estimating the DWA spurs using stimulus and DNL error shape:
[n_hat,ck,E] = dwa_spur_est(Bits,BPmode,Estd,shape,k,Bf,data,IDWAmode,1);
%n_hat is the estimate

if random_dir
    seqx = sign(randn(1,DL));
    seqx(seqx==0)=1;
else
    seqx=ones(1,DL);
end

if tf_mode
    seq = seqx;  % ...for plain DWA
else
                    % ...for BNF-dwa (see the ESL block)
    seq = [seqx -seqx]; 
end

if DACmode
    if IDWAmode
        sim IDWAMOD_DAC
    else
        sim DWAMOD_DAC
    end
else
    if IDWAmode
        sim IDWAMOD_ADC
    else
        sim DWAMOD_ADC
    end
end

%---- PLOT
if not(exist('win','var')) || (exist('win','var') && length(win)~=NT)
    %win=kaiser(NT,20);
    win=hann(NT);
end

if BPmode
    NS=pow2(log2(NT)-10);
else
    NS=pow2(log2(NT)-9);
end
wins = hann(NT/NS);

HdB=20*log10(abs(fft(win.*simout,NT)/(sum(win)*0.5)));
MdB=20*log10(abs(fft(win.*mmnoise,NT)/(sum(win)*0.5)));
MestdB=20*log10(abs(fft(win.*n_hat.',NT)/(sum(win)*0.5)));

Hid_dB=20*log10(abs(fft(win.*(simout+mmnoise),NT)/(sum(win)*0.5)));
fvek=linspace(fs/NT,fs/2,NT/2);

figure
subplot(321)
plot(fvek,HdB(1:(NT*0.5)),'b')
title('output signal'),axis tight
xlabel('frequency'),ylabel('dB')
hold on
YLim=get(gca,'YLim');
if BPmode
    plot((fs/4-fb/2)*[1 1],YLim,'k--',(fs/4+fb/2)*[1 1],YLim,'k--')
    text(fs/4+fb/2,YLim(2)-10,'\leftarrowsignal band')
else
    plot([fb fb],YLim,'k--')
    text(fb,YLim(2)-10,'\leftarrowsignal band')
end
hold off

subplot(325)
plotflag = 1;
%switch modelname
%    case {'DSDAC_none','DSADC_none'}
%        plot(fvek,MdB(1:(NT*0.5)),'r')
%        title('mismatch noise'),axis tight
%        xlabel('frequency'),ylabel('dB')
%        
%    otherwise
        plot(fvek,MestdB(1:(NT*0.5)),'c',fvek,MdB(1:(NT*0.5)),'r')
        title('mismatch spur estimate & true mismatch noise'),axis tight
        xlabel('frequency'),ylabel('dB')
%end
hold on        
if BPmode
    plot((fs/4-fb/2)*[1 1],YLim,'k--',(fs/4+fb/2)*[1 1],YLim,'k--')
    text(fs/4+fb/2,YLim(2)-10,'\leftarrowsignal band')
else
    plot([fb fb],YLim,'k--')
    text(fb,YLim(2)-10,'\leftarrowsignal band')
end
hold off
set(gca,'YLim',YLim)

subplot(323)
plot(fvek,Hid_dB(1:(NT*0.5)),'k')
title('ideal output'),axis tight
xlabel('frequency'),ylabel('dB')
hold on
%YLim=get(gca,'YLim');
if BPmode
    plot((fs/4-fb/2)*[1 1],YLim,'k--',(fs/4+fb/2)*[1 1],YLim,'k--')
    text(fs/4+fb/2,YLim(2)-10,'\leftarrowsignal band')
else
    plot([fb fb],YLim,'k--')
    text(fb,YLim(2)-10,'\leftarrowsignal band')
end
hold off
set(gca,'YLim',YLim)

%Subplots below use Modulo Time Plot, see File ID: #22907
t=linspace(0,Tfinal,NT).';

subplot(322)
tt=mod(t,1/fsig);
[tt,index]=sort(tt);
yy=simout(index);
n_poly = 5;
warning off %#ok<WNOFF>
p=polyfit(tt,yy,n_poly);
warning on %#ok<WNON>
yest=polyval(p,tt);
plot(tt,yy,'b.',tt,yest,'--k','MarkerSize',4)
title('modulo time plot, output signal')
axis tight,legend('data','mean',1),
xlabel(['time, 1 period is ', num2str(1/fsig,2),'s. t_s = ',...
    num2str(ts,2),'s'])

subplot(324)
qq = yy-yest;
plot(tt,qq,'m.','MarkerSize',4)
title('modulo time plot, qv. noise'), axis tight,xlabel('time, 1 period')
subplot(326)
plot(tt,mmnoise(index),'r.','MarkerSize',4)
axis tight, title('modulo time plot, mismatch noise')
xlabel('time, 1 period')

figure
subplot(211)
spectrogram(mmnoise,wins,length(wins)*0.5,[],1)
title('true mismatch noise spectrogram')
subplot(212)
spectrogram(n_hat,wins,length(wins)*0.5,[],1)
title('mismatch noise spurious noise estimate')

Contact us