12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273 |
- function H = spm_dwtmtx(N,K,T)
- % Create basis functions for Discrete (Haar) Wavelet Transform
- % FORMAT H = spm_dwtmtx(N,K,T)
- %
- % N - dimension
- % K - order: number of basis functions = N/K
- %
- % T - option flag for thinning eccentric wavelets [default: false]
- %__________________________________________________________________________
- %
- % spm_dwtmtx creates a matrix for the first few basis functions of a one
- % dimensional Haar Discrete Wavelet transform.
- %__________________________________________________________________________
- % Copyright (C) 2011-2015 Wellcome Trust Centre for Neuroimaging
-
- % Karl Friston
- % $Id: spm_dwtmtx.m 6416 2015-04-21 15:34:10Z guillaume $
-
-
- % Create basis set
- %==========================================================================
- try, T; catch, T = false; end
-
- K = max(K,1);
- H = ones(N,1);
- I = 1;
- while ~isempty(I)
-
- % Find non-zero elements and create new bases
- %----------------------------------------------------------------------
- j = find(H(:,I(1)) > 0);
- k = find(H(:,I(1)) < 0);
- nj = length(j);
- nk = length(k);
- hj = sparse(1:fix(nj/2),1,2,nj,1) - 1;
- hk = sparse(1:fix(nk/2),1,2,nk,1) - 1;
-
- % Add new columns if there are enough elements
- %----------------------------------------------------------------------
- if nj > K
- H(j,end + 1) = hj;
- I(end + 1) = size(H,2);
- end
- if nk > K
- H(k,end + 1) = hk;
- I(end + 1) = size(H,2);
- end
-
- % This column has now been expanded
- %----------------------------------------------------------------------
- I(1) = [];
-
- end
-
- % Thin eccentric basis functions
- %==========================================================================
- if T
-
- % indices of column (bases) to retain
- %----------------------------------------------------------------------
- j = logical(H(:,1));
- for i = 1:length(j);
- k = find(H(:,i));
- l = length(k);
- if any(k < (N/2 - 2*l)) || any(k > (N/2 + 2*l))
- j(i) = 0;
- end
- end
-
- % thin
- %----------------------------------------------------------------------
- H = H(:,j);
- end
|