function c = multSIDt(sid, y)
% multSIDt        Multiplication by transposed Shift Invariant Dictionary: c = D'*y
%                 
% Except for very small SID, fft is used to do the calculations.
%
% use:   c = multSIDt(sid, y);
% -------------------------------------------------------------------------
%   c     a vector Kx1 , or matrix KxL, 
%         innerproducts of dictionary atoms and signal y
%   sid   The Shift Invariant Dictionary, a struct as in initSID
%   y     a vector Nx1 , or matrix NxL
% -------------------------------------------------------------------------

%----------------------------------------------------------------------
% 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  29.10.2015  KS: made file 
%----------------------------------------------------------------------

if ~isstruct(sid) || ~isfield(sid,'Type') || ~strcmp(sid.Type, 'SID')
    disp([mfilename,': First argument not a SID structure.']);
    c = [];
    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);
[N,L] = size(y);

if N < Qmax
    disp([mfilename,': signal length smaller than longest atom.']);
    c = [];
    return
end

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


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

if N*K < 1e4    
    % I found that using fft is almost always faster!
    % this may be a good test to wether fft gives correct results or not
    D = makeSIDmatrix(sid, N);
    c = D'*y;
else      % use fft 
    c = zeros(K,L);
    NN = 2^ceil(log2(N + Qmax - 1));
    ffty = fft(y,NN);
    for m = 1:M
        atom = sid.freeVar(sid.i1(m):sid.i2(m));
        fftatom = fft(flipud(atom), NN);
        cm = ifft(ffty .* repmat(fftatom,1,L));
        % fprintf('Atom %i has length %i, NN = %i, size(cm) = %i x % i \n', ...
        %      m, numel(atom), NN, size(cm,1), size(cm,2) );
        % fprintf('Update c(%i:%i) using cm(%i:%i:%i) \n', ...
        %     k1(m), k2(m), Qm(m),Pm(m),(Qm(m)+(k2(m)-k1(m))*Pm(m)));
        %
        c(k1(m):k2(m),:) = cm(Qm(m):Pm(m):(Qm(m)+(k2(m)-k1(m))*Pm(m)),:);
    end
end
    
end






