create_DWT

 Creates a DWT representation matrix without signal extension

0001 % Creates a DWT representation matrix without signal extension
0002 
0003 function Psi = create_DWT(in)
0004 
0005 % A wavelet basis function using orthogonal wavelets
0006 %
0007 % inputs:
0008 %   J -- finest scale for wavelets
0009 %   wType -- type of wavelets
0010 %   sym -- type of extension
0011 %   N -- length of signal interval
0012 
0013 N = in.N;  
0014 wType = in.wType; 
0015 J = in.J;
0016 sym = in.sym;
0017 
0018 if strcmpi(wType(1:4),'daub')
0019     switch wType
0020         case 'daub79'
0021             [h0 h1 g0 g1] = daub79; sym = 1;
0022         case 'daub1018'
0023             [h0 h1 g0 g1] = daub1018; sym = 2;
0024             %         case 'daub4'
0025             %             filter_length = 4;
0026             %             [h0 h1 g0 g1] = daub4;
0027             %         case 'daub8'
0028             %             filter_length = 8;
0029             %             [h0 h1 g0 g1] = daub8;
0030         otherwise
0031             filter_length = str2num(wType(5:end));
0032             [h0 h1 g0 g1] = dauborth(filter_length);
0033     end
0034 else
0035     error('Not implemented yet');
0036 end
0037 
0038 % switch lower(wType)
0039 %     case 'daub2'
0040 %         [h0 h1 g0 g1] = dauborth(2);
0041 %         filter_length = 2;
0042 %     case 'daub4'
0043 %         [h0 h1 g0 g1] = dauborth(4);
0044 %         filter_length = 4;
0045 %     case 'daub8'
0046 %         [h0 h1 g0 g1] = dauborth(8);
0047 %         filter_length = 8;
0048 %     case 'daub12'
0049 %         [h0 h1 g0 g1] = dauborth(12);
0050 %         filter_length = 12;
0051 %     otherwise
0052 %         disp('NAO');
0053 % end
0054 
0055 % compute the lengths of functions at every scale
0056 if sym == 3
0057     w_len = filter_length;
0058     % for j = 1:J-1; w_len(end+1) = w_len(end)+2^j*(filter_length-1); end
0059     for j = 1:J-1; w_len(end+1) = 2*w_len(end)-1+(filter_length-1); end
0060     w_len(end+1) = w_len(end); 
0061 else
0062     w_len = N;
0063 end
0064 %
0065 if w_len(end) > N && sym == 3;
0066     error('reduce the scale J or increase length N - length of scaling function is larger than N');
0067 end
0068 
0069 % compute synthesis matrix...
0070 switch sym
0071     case 3
0072         iW_h = @(x) idwtmult1_conv(x,g0,g1,J);
0073         % iW_h = @(x) dwtmult1_conv(x,h0,h1,J);
0074     otherwise
0075         iW_h = @(x) ifwt(x,g0,g1,J,sym);
0076 end
0077 Psi = [];
0078 for ii = 1:N;   
0079     Psi(:,end+1) = iW_h(circshift([1; zeros(N-1,1)],ii-1));    
0080 end
0081 
0082 % %-----------------% %
0083 % % debug
0084 % %-----------------% %
0085 % PSI = [[Psi; zeros(N)] [zeros(N); Psi]];
0086 % vec = @(z) z(:);
0087 % max(abs(vec(PSI'*PSI-eye(2*N))))
0088 % %%
0089 % in = []; in.J = 4; in.N = 256; in.wType = 'daub10'; in.sym = 3;
0090 % Psi = create_DWT(in);