dwa_noise_model(N,BPmode,Estd,shape,k,u,plotflag)

function [mm_hat,E] = dwa_noise_model(N,BPmode,Estd,shape,k,u,plotflag)
% DWA DAC mismatch noise model 
% u      DSM stimulus normalized between [0,N] 
%           (just like in simulateESL in DS toolbox)
% N         The number of unit DACs
%           N can be incremented to model Incremental DWA: different scaling
%           for u
% Bpmode    1, if bandpass DWA, else 0
%           note that u should be bandpass as well for BPmode = 1
% Estd      stdev of DNL error, the mean error will be 0
%           OR
%           for custom shape 'cust', Estd is the DNL error vector
% shape     'poly' for polynomial shape
%           'sin' for sinusoida shape
%           'rand' for random shape
%           'cust' for custom shape
% k         polynomial order or
%           periods in sinusoid
%           for 'rand' and 'cust' shape, k is negligible
% plotflag  if exists and is 1, plot (else, don't)
% mm_hat    shaped mismatch noise estimate
% E         The DNL error vector
% 
% examples:
%
% polynomial DNL error vector shape
% [mm_hat,ck,E,Af,Bf] = dwa_spurs(N,BPmode,Estd,'poly',k,u,plotflag)
% [mm_hat,ck,E,Af,Bf] = dwa_spurs(7,0,0.01,'poly',3,3*ones(1,256),1);
%
% sinusoidal DNL
% [mm_hat,ck,E,Af,Bf] = dwa_spurs(N,BPmode,Estd,'sine',k,u,plotflag)
% [mm_hat,ck,E,Af,Bf] = dwa_spurs(7,0,0.01,'sine',1,3*ones(1,256),1);
%
% custom DNL
% [mm_hat,ck,E,Af,Bf] = dwa_spurs(N,BPmode,Estd,E,0,u,plotflag)
% [mm_hat,ck,E,Af,Bf] = dwa_spurs(7,0,0.01,[6,1,3,6,6,7,5],0,3*ones(1,256),1);
%
% random DNL (normal distribution)
% [mm_hat,ck,E,Af,Bf] = dwa_spurs(N,BPmode,Estd,E,'rand',u,plotflag)
% [mm_hat,ck,E,Af,Bf] = dwa_spurs(7,0,0.01,'rand',0,3*ones(1,256),1);

if nargin == 6
    plotflag = 0;
end

NT = length(u);

if isnumeric(shape)
    E = shape;
    shape = 'cust';
end

switch shape
    case 'poly'
        polyflag = 1;  % if 1, polynomial DNL
    otherwise
        polyflag = 0;
end

if N < 4
    Nspurs = N-1; 
else
    Nspurs = floor(N/2);
end

% select spur(s) at frequency k*x/N. 
% if polyflag, k = polynomial order
% if k = 0, random
switch shape 
    case 'sine'
        if k > Nspurs
            fprintf('Max value is %i for LEVS = %i',Nspurs,round(N/2))
            fprintf('\nQuitting')
            return
        end
    case 'poly'
        if k == 0
            fprintf('Min value for k is 1')
            fprintf('\nQuitting')
        end
    case 'rand'
        E = randn(1,N);
        shape = 'cust';
end



fs=1;

if BPmode
    fc = fs/4;
    Af = [1 ];
    Bf = [1 0 1];
else
    fc = 0;
    Af = 1;
    Bf = [1 -1];
end

%  DNL begin
switch shape
    case 'cust'
        %E = dsm_dnl(N,mean(u));
        if length(E)~=N
            fprintf('The length of DNL error is %i. Should be %i',length(E),N)
            fprintf('\nQuitting\n')
            return
        else
            fprintf('Custom shape used and scaled in stdev of %g\n',Estd)
            
        end
    otherwise
        E = zeros(1,N);
        if any(k==0)
            E = randn(size(E));
        else
            if polyflag
                E=real((linspace(-1,1,N)).^k);
    	
            else
                E  = E+ sin(k*(0:(N-1))/(N)*2*pi);
            end
        end
end
E = E-mean(E);
E = Estd*E/std(E);

EX = fft(E.');
EXmag = abs(EX);
ck = EXmag(2:(Nspurs+1));
EXang = angle(EX)/pi*180;
% DNL end

%y = zeros(Nspurs,NT);
mm_hat = zeros(1,NT);
beta1 = 1/pi;

p=[-1.958948344182117   5.044492106506190  -7.055295364055410   8.334364561858429  -5.365817455286532  -0.064285846600190   2.257543740766184];
qg = polyval(p,abs(0.5*u/N));

if BPmode
    phi = (1-qg)*2*pi*2;
else
    phi = (1-qg)*2*pi;
end

switch shape
    case 'sine'
        fd=k/N;
        mm_hat = beta1*ck(k)/k*cos(fc*2*pi + 2*pi*fd*(filter(Af,Bf,u)));
    otherwise
        for k = 1:Nspurs
            fd=k/N;  
            mm_hat = mm_hat + beta1*ck(k)/k*cos(fc*2*pi + 2*pi*fd*(filter(Af,Bf,u))+phi);
        end
end

%the shaping
if BPmode
    mm_hat = filter(Bf,Af,mm_hat);
else
    mm_hat = filter(Bf,Af,mm_hat);
end

if plotflag
    f_E = 0:1/N:1;
    f_E(end)=[];

    xtl = cell(1,N);
    ktl = cell(1,N);
    for ind = 1:length(f_E)
        temp = rats(f_E(ind));
        temp(findstr(temp,' '))=[];
        xtl(ind) = {char(temp)};
        ktl(ind) = {['k = ',int2str(ind)]};
    end
    %figure
    subplot(311)
    plot(E,'-o')
    ylabel('\epsilon')
    title('DNL error vector')
    subplot(3,2,[3 5])
    stem(f_E,EXmag)
    set(gca,'XTick',f_E)
    set(gca,'XLim',[0,1])
    set(gca,'XTickLabel',xtl)
    ylabel('|fft( \epsilon )|')
    xlabel('freq/fs')

    subplot(3,2,[4 6])
    stem(f_E,EXang)
    set(gca,'XTick',f_E)
    set(gca,'XLim',[0,1])
    set(gca,'XTickLabel',xtl)
    set(gca,'YTick',-180:45:180)
    ylabel('angle(fft( \epsilon ))')
    xlabel('freq/fs')
end