function sid = initSID(Qm, Pm, y)
% initSID         Initialize Shift Invariant Dictionary given atom lengths and shifts.
%
% If no example signal is given as input, atoms are initialized by fir1.
% In this case the longest atoms are put first, and the atoms are also 
% sorted in frequency from low-pass to high-pass. Initial values are 
% from Matlab fir1 function for filter design.
%
% If an example signal is given (three input arguments) initial values are
% made from this signal.
%
% use:   sid = initSID(Qm);
%        sid = initSID(Qm, Pm);
%        sid = initSID(Qm, Pm, y);
% -------------------------------------------------------------------------
%   Qm    The length of each atom, size Mx1
%   Pm    The shift to use for each atom, size Mx1
%         if not given, all shifts are assumed to be 1
%   y     An example signal, size Nx1, N should be in range 1000-20000
%         or the values to use, size Qx1,  Q = sum(Qm)
%         (N > Q, only works if dictionary learning is available)
%   sid   The Shift Invariant Dictionary, a struct with fields
%         Type     : always  'SID'
%         Text     : always  'Shift Invariant Dictionary'
%         M        : number of atoms in dictionary
%         Q        : total number of free variables, i.e. sum of Qm
%         Qm       : length of each atom, size Mx1
%         Pm       : shift for each atom, size Mx1
%         freeVar  : the free variables, a vector of size Qx1
%                    (name: q, qvec, x, or xvec instead???)
%         i1,i2    : perhaps also these, atom m is 
%                    sid.freeVar(sid.i1(m):sid.i2(m))
%         Other fields may be included when needed, ex:
%         AtAc     : compact (QxM) representation of A'*A, for learning
%         Atb      : representation of A'*b for dictionary learning
%         Info     : text telling how dictionary was learned/used
% -------------------------------------------------------------------------
% ex:   Qm = [64,32,16,8,8]; Pm = [8,4,2,1,1]; sid = initSID(Qm,Pm);
%       Qm = [35,25,20,10,10];  Pm = [ 5, 2, 2, 1, 1]; sid = initSID(Qm,Pm);
%       sid = initSID(Qm, Pm, randn(sum(Qm),1));

%----------------------------------------------------------------------
% 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 
% Ver. 1.1  17.11.2015  added field 'Info'
% Ver. 1.2  18.11.2015  added third input argument
%           20.11.2015  and made it more flexible
%----------------------------------------------------------------------


M = numel(Qm);
Q = sum(Qm);
if (nargin < 2) || (numel(Pm) ~= M)
    Pm = ones(M,1);
end

i2 = cumsum(Qm(:));
i1 = [1; i2(1:(M-1))+1];

Qm = Qm(:);
Pm = Pm(:);
sid = struct( 'Type', 'SID', ...
    'Text', 'Shift Invariant Dictionary', ...
    'Info', ['Initialized by ',mfilename,' ',datestr(now),'.'], ...
    'M', M, ...      % number of atoms
    'Q', Q, ...      % total number of free variables
    'Qm',Qm, ...     % length of each atom
    'Pm',Pm, ...     % shift for each atom
    'freeVar',zeros(Q,1), ... % the free variables
    'i1',i1, ...
    'i2',i2  );

if (nargin <= 2)   % initialize using fir1
    %
    [Qm, I] = sort(Qm,'descend');
    Pm = reshape(max(Pm(I),1),M,1);
    i2 = cumsum(Qm(:));
    i1 = [1; i2(1:(M-1))+1];
    sid.Qm = Qm;
    sid.Pm = Pm;
    sid.i1 = i1;
    sid.i2 = i2;
    
    if M == 1
        sid.freeVar = fir1(Q-1, 0.5)';
    elseif M < 6    % few atoms use dyadic split of frequencies
        m = M;  % last atom is high pass
        if rem(Qm(m),2)  % odd
            sid.freeVar(i1(m):i2(m)) = fir1(Qm(m)-1, 0.5, 'high')';
        else  % even
            sid.freeVar(i1(m):i2(m)) = [fir1(Qm(m)-2, 0.5, 'high'),0]';
        end
        f = [0.25, 0.5];
        for m = (M-1):(-1):2
            % middle atoms are band pass
            sid.freeVar(i1(m):i2(m)) = fir1(Qm(m)-1, f)';
            f = f/2;
        end
        m = 1;   % first atom is low-pass
        sid.freeVar(i1(m):i2(m)) = fir1(Qm(m)-1, f(2))';
        %
    else  % more atoms use uniform split of frequencies
        m = M;
        f = [1-1/M, 1];
        if rem(Qm(m),2)  % odd
            sid.freeVar(i1(m):i2(m)) = fir1(Qm(m)-1, f(1), 'high')';
        else  % even
            sid.freeVar(i1(m):i2(m)) = [fir1(Qm(m)-2, f(1), 'high'),0]';
        end
        for m = (M-1):(-1):2
            % middle atoms are band pass
            f = f - 1/M;
            sid.freeVar(i1(m):i2(m)) = fir1(Qm(m)-1, f)';
        end
        m = 1;   % first atom is low-pass
        sid.freeVar(i1(m):i2(m)) = fir1(Qm(m)-1, f(1))';
    end
    sid.Info = ['Initialized by ',mfilename,' ',datestr(now),' using FIR filters.'];
    
else  % three or more input arguments, initialize from input signal y
    %
    y = y(:);
    N = numel(y);
    if N == 1
        sid.freeVar = y(1)*ones(Q,1);
        sid.Info = ['Initialized by ',mfilename,' ',datestr(now),' setting all entries to ',num2str(y(1)),'.'];
    elseif (N < Q)
        sid.freeVar = [y; ones(Q-N,1)];
        sid.Info = ['Initialized by ',mfilename,' ',datestr(now),' using supplied values and ones.'];
    elseif (N == Q)
        sid.freeVar = y(:);
        sid.Info = ['Initialized by ',mfilename,' ',datestr(now),' using supplied values.'];
    else  % N > Q
        %
        r = y;
        h = [1,3,3,1];   % use a short low pass filter
        for i = 1:M    % find atoms one by one
            sidi = struct( 'Type', 'SID', ...
                'Text', 'Shift Invariant Dictionary', ...
                'Info', 'temporary SID.', ...
                'M', i, ...
                'Q', sum(Qm(1:i)), ...
                'Qm',Qm(1:i), ...
                'Pm',ones(i,1), ...
                'freeVar',zeros(sum(Qm(1:i)),1), ...
                'i1',i1(1:i), ...
                'i2',i2(1:i)  );
            if (i > 1) % keep previous free variables
                sidi.freeVar(1:i2(i-1)) = freeVar;  % from last
            end
            aL = Qm(i);  % atom length
            [temp,k1] = max(filter(h,1,abs(r)));
            k1 = min(max(1,k1-floor(aL/2)),N-aL);
            k2 = k1 + aL - 1;
            atom = r(k1:k2);
            na = norm(atom);
            sidi.freeVar(i1(i):i2(i)) = atom/na;
            if ((exist('learnSID','file') == 2) && (exist('saSIDbmp_mex','file') == 3))
                p1 =  struct('tSNR',40, 'tsf',i/(10*M), 'lambda0',0.5, 'lambda1',0.5, 'noIt',5);
                sidi = learnSID(y, sidi, p1, 0);
                s = floor(p1.tsf*N);
                w = sparse(saSIDbmp_mex(y, sidi.freeVar, sidi.Qm, sidi.Pm, s, 0));
            else
                s = floor((N*i)/(10*M));
                w = saSIDbmp2(y, sidi, s);
            end
            yr = multSID(sidi, w, N);
            r = y - yr;
            freeVar = sidi.freeVar;
        end
        sid.freeVar = freeVar;
        sid.Info = ['Initialized by ',mfilename,' ',datestr(now),' using atoms learned from example signal.'];
    end
    % 
end

% normalize
for m = 1:M
    a = sid.freeVar(i1(m):i2(m));
    sid.freeVar(i1(m):i2(m)) = a/norm(a);
end

end

