function y = multSID(sid, w, N)
% multSID         Multiplication for Shift Invariant Dictionary: y = D*w
%
% use:   y = multSID(sid, w, N);
% -------------------------------------------------------------------------
%   y     a vector Nx1 , or matrix NxL
%   sid   The Shift Invariant Dictionary, a struct as in initSID
%   w     a sparse vector Kx1 (or matrix, KxL)
%   N     length of y
% -------------------------------------------------------------------------
% ex:   Qm = [8,8,16,32,64]; Pm = [1,1,2,4,8]; sid = initSID(Qm,Pm);

%----------------------------------------------------------------------
% Made by Karl Skretting, University of Stavanger (UiS), 
% Department of Electrical Engineering and Computer Science (IDE)
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  28.10.2015  KS: made file 
%----------------------------------------------------------------------

if ~isstruct(sid) || ~isfield(sid,'Type') || ~strcmp(sid.Type, 'SID')
    disp([mfilename,': First argument not a SID structure.']);
    y = [];
    return
end
if nargin < 3
    disp([mfilename,': use N (length of output) as third argument.']);
    y = [];
    return
end

M = sid.M;    %  number of atoms in dictionary
% Q = sid.Q;    %  total number of free variables, i.e. sum of Qm
Qm = sid.Qm;  %  length of each atom, size Mx1
Pm = sid.Pm;  %  shift for each atom, size Mx1
% Qmax = max(Qm);

Km = floor( (N+Pm-Qm)./Pm );
K = sum(Km);
k2 = cumsum(Km);
k1 = [1; k2(1:(M-1))+1];

% check length of K
if K ~= size(w,1)
    % w is size KxL, and this size should match the K value calculated from
    % sid and N, for each N value the correct K value can be found by
    %   fK = @(N) sum(floor((N+Pm-Qm)./Pm));         
    % or as used if we want it to be called ny a row vector of N-values
    %   fK = @(N) sum(floor( (repmat(reshape(N,1,numel(N)),numel(Pm),1) + ...
    %          repmat(Pm(:),1,numel(N))-repmat(Qm(:),1,numel(N)))./ ...
    %          repmat(Pm(:),1,numel(N)) ));   
    % The inverse of this is not as straighforwared as wanted, could be
    %   fN = @(K) ceil((K + sum(Qm./Pm))/sum(1./Pm));     
    % A fixed rule for how to do the calculation of N, given K could be
    % made. Observe that the increase in fK is periodic with periode as
    % least common multiple of numbers in Pm. 
    % Here we do it the 'simple' way
    disp([mfilename,': value of N and size of w do not match.']);
    y = [];
    return
end

% many methods can be used to calculate y
% depending on sparseness and size of w it may vary which method is best

sf = nnz(w)/numel(w);   % NOTE: this is not nnz(w)/numel(y)
if ((N*K < 1e7) && (sf > 0.1)) || ((N*K < 1e5) && (sf > 0.02))  
    % this is often faster
    D = makeSIDmatrix(sid, N);
    y = D*w;
else      % this is good for quite sparse w   nnz(w)/numel(w) < 0.15
    L = size(w,2);
    y = zeros(N,L);
    for i = 1:L
        coeff = nonzeros(w(:,i));
        idx = find(w(:,i));
        m = 1;
        atom = sid.freeVar(sid.i1(m):sid.i2(m));
        for j = 1:numel(idx)
            k = idx(j);
            while k > k2(m)
                m = m + 1;
                atom = sid.freeVar(sid.i1(m):sid.i2(m));
            end
            k = k - k1(m);   % k is column in Dm part of D;  0,1,...Qm(m)-1
            n1 = k*Pm(m)+1;  % start for n
            n = n1:(n1+Qm(m)-1);
            y(n,i) = y(n,i) + coeff(j)*atom;
        end
    end
end

if nnz(y) > 0.5*numel(y)
    y = full(y);
else
    y = sparse(y);
end

end






