% Create a synthesis matrix containing mexican hat wavelets 
%
% waveletmat.m
%
% inputs: N ......... length of signal
%         dilparam .. scale 'offset' parameter or vector of scales
%         shift ..... shift parameter
%         w_type .... wavelettype: 
%                     0 .. first scaling then shifting
%                          g_k,l = T_ls D_kd g
%                          corresponds to 'Pyramid algorithm'
%                     1 .. first shifting then scaling
%                          g_k,l = D_kd T_ls  g
%                          corresponds to 'Algorithm à trous'
%
% notes: The resulting matrix corresponds to the synthesis matrix of the
%        wavelet sequence using the standard mexican hat functions as 
%        mother wavelet. In the rows of this matrix the single elements 
%        can be found such that:
%        D_(ka) T_(lb) g = WM[l+k*log(N)/log(dilparam)] 
%           (for w_type = 1 and dilparam a number)
%
% references:
%        [Mallat] S. Mallat, A Wavelet Tour of Signal Processing
%        [XXL] P. Balazs, Hilbert-Schmidt Operators and Frames - 
%              Classification, Best Approximation by Multipliers and 
%              Algorithms!
%
% author: XXL .. peter.balazs@oeaw.ac.at
%
% Thanks to Laurent Jacques
%
% See also: TestBestApprMult.m
%
% Copyright :  (c) NUHAG, Faculty of Math., University of Vienna, AUSTRIA
%              http://nuhag.mat.univie.ac.at/
%              (c) Acoustics Research Institute, Austrian Academy of
%              Science
%              http://www.kfs.oeaw.ac.at
function [WM xpo] = waveletmat(N,dilparam,shift,w_type);

if nargin < 1
    N = 256;
end
if nargin < 2
    dilparam = 2;
else
    [Mdil Ndil] = size(dilparam);
    if Mdil ~= 1 && Ndil ~= 1
        error('Parameter dilparam has to be a number or a vector');
        return;
    elseif Mdil ~= 1 || Ndil ~=1
        a = dilparam;
        M = length(a);
    else
        M = ceil(log(N)/log(dilparam))-1; % number of scales
        a = [dilparam.^(-M:-1) 1 dilparam.^(1:M)];
    end
end
if nargin < 3
    shift = 1;
end
if nargin < 4
    w_type = 0;
end

t = 0:(N-1);

WM = [];
xpo = [];
for k = 1:length(a)
   mother = (1/a(k))^.5 * mexicanhat((t - N/2)/a(k));
   if w_type == 0
       s_shift = shift;
   else
       s_shift = max(ceil(shift*a(k)),shift);
   end
   sss = 1:s_shift:N;
   aaa = zeros(1,N/s_shift);
   aaa(1:(N/s_shift)) = a(k);
   xpo = [xpo [sss ; aaa ]];
   WM = [WM; circul(mother,s_shift)];
end