dwa_spur_est(Bits,BPmode,Estd,shape,k,Bf,data,IDWAmode,plotflag)

function [n_hat,ck,E,N] = dwa_spur_est(Bits,BPmode,Estd,shape,k,Bf,data,IDWAmode,plotflag)

if nargin < 7
    fprintf('No demo available\n')
    return
elseif nargin == 7
    IDWAmode = 0;
    plotflag = 1;
elseif nargin == 8
    plotflag = 1;
end

NT = length(data);

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

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

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

Nspurs_default = floor(N/2);
Nspurs = Nspurs_default;

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



fs=1;
Af = 1;

if BPmode
    fc = fs/4;

    if IDWAmode
        fprintf('\nNote: IDWA is not advantageous in bandpass mode.\n')
    %    IDWAmode = 0;
    end
else
    fc = 0;
end



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

%  DNL begin
switch shape
    case 'cust'
        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 with 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);   %EX.*conj(EX)/N;
ck = EXmag(2:(Nspurs+1));
EXang = angle(EX)/pi*180;
% DNL end

y = zeros(Nspurs,NT);

% constant amplitudes
x = WLoffset+WLscale*data;

for k = 1:Nspurs
    fd=k/N;
    %y(k,:) = EXmag(k+1)*cos(fc*2*pi + 2*pi*fd*mod(filter(Af,Bf,x),N))/k;
    y(k,:) = 0.5*ck(k)/k*cos(fc*2*pi + 2*pi*fd*mod(filter(Af,Bf,x),N));    
end
n_hat = sum(y);

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')

%     fvec = linspace(0,0.5,NT/2);
%     klegend = cell(1,Nspurs);
%     for k = 1:Nspurs
%         klegend(k) = {['k = ' int2str(k)]};
%     end
% 
%     figure
%     Y = zeros(Nspurs,NT);
%     for k = 1:Nspurs
%         Y(k,:) = abs(fft(y(k,:)))/NT;
%     end
%     subplot(211)
%     plot(fvec,(Y(:,1:(NT/2))))
%     title('DWA spurs magnitude spectrum')
%     legend(klegend,2)
%     set(gca,'XLim',[0,0.5])
%     
%     
%     if BPmode
%         NS=pow2(log2(NT)-10);
%     else
%         NS=pow2(log2(NT)-9);
%     end
%     wins = kaiser(NT/NS,20);
%     subplot(212)
%     spectrogram(n_hat,wins,length(wins)*0.5,[],1)
%     title('DWA spurs spectrogram')
%     set(gca,'XLim',[0,0.5])
end