## [77a22a]: wavelets / fwt.m  Maximize  Restore  History

### 158 lines (143 with data), 6.1 kB

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155``` ```function c = fwt(f,h,J,varargin) %FWT Fast Wavelet Transform % Usage: c = fwt(f,h,J); % c = fwt(f,h,J,...); % % Input parameters: % f : Input data. % h : Analysis Wavelet Filterbank. % J : Number of filterbank iterations. % a : Explicitly defined subsampling factors. % % Output parameters: % c : Coefficients stored in a cell-array. % % `c=fwt(f,h,J)` returns wavelet coefficients *c* of the input signal *f* % using *J* iterations of the basic wavelet filterbank defined by *h*. % The fast wavelet transform algorithm (or Mallat's algorithm) is employed. % If *f* is a matrix, the transformation is applied to each of *W* columns. % % The coefficents *c* are Discrete Wavelet transform of the input signal *f*, % if *h* defines two-channel wavelet filterbank. The function can apply % the Mallat's algorithm using basic filter banks with any number of the % channels. In such case, the transform have different name. % % The basic analysis wavelet filterbank \$h\$ can be passed in several formats. % The formats are the same as for the |fwtinit|_ function. The simplest % is passing a cell array, whose first element is the name of the function % defining the basic wavelet filters (`wfilt_` prefix) and the other elements are % the parameters of the function. e.g. `{'db',10}` calls `wfilt_db(10)` internally. % The second possible format of \$h\$ is to pass cell array of one dimensional % numerical vectors directly defining the wavelet filter impulse responses. % In this case, outputs of the filters are subsampled by a factor equal % to the number of the filters. For creating completely custom filterbanks use the % |fwtinit|_ function. % The third option is to pass structure obtained from the |fwtinit|_ % function. % % The coefficients in the cell array \$c\{jj\}\$ for \$jj=1,\ldots,J\cdot(N-1)+1\$, % where \$N\$ is number of filters in the basic wavelet filterbank, are ordered % with inceasing central frequency of the corresponding effective filter frequency % responses or equivalently with decreasing wavelet scale. % % If the input *f* is matrix with *W* columns, each element of the cell % array \$c\{jj\}\$ is a matrix with *W* columns with coefficients belonging % to the appropriate input channel. % % The following flag groups are supported: % % 'per','zpd','sym','symw','asym','asymw','ppd','sp0' % Type of the boundary handling. % % % Boundary handling: % ------------------ % % The default periodic extension considers the input signal as it was a one % period of some infinite periodic signal as is natural for the transforms % based on the FFT. The resulting wavelet representation is % non-expansive, that is if the input signal length is a multiple of a % \$J\$-th power of the subsampling factor, the total number of coefficients % is equal to the input signal length. % % If the input signal length is not a multiple of a \$J\$-th power of the % subsampling factor, the processed signal is padded internally by % repeating the last sample at each step of the transform to the next % multiple of the subsampling factor rather than doing the prior explicit % padding. In addition, the periodic extension restrict the input signal % length to be greater than a certain length. % % `fwt(f,h,J,ext)` with `ext` other than `'per'` computes a slightly % redundant wavelet representation of the input signal *f* with the chosen % boundary extension *ext*. % % The default periodic extension can result in "false" high wavelet % coefficients near the boundaries due to the possible discontinuity % introduced by the periodic extension. The custom extension can diminish % this phenomenon. The extensions are done at each level of the transform % internally rather than doing the prior explicit padding. The supported % possibilities are: % % * `'zpd'` - zeros are considered outside of the signal (coefficient) support. % % * `'sym'` - half-point symmetric extension. % % * `'symw'` - whole-point symmetric extension % % * `'asym'` - half-point antisymmetric extension % % * `'asymw'` - whole point antisymmetric extension % % * `'ppd'` - periodic padding, same as `'ppd'` but the result is expansive representation % % * `'sp0'` - repeating boundary sample % % Note that the same flag have to be used in the call of the inverse transform % function `ifwt`. % % Examples: % --------- % % A simple example of calling the |fwt|_ function::: % % f = gspi; % J = 10; % c = fwt(f,{'db',8},J); % plotfwt(c); % % See also: ifwt, fwtinit % % References: ma98 if nargin<3 error('%s: Too few input parameters.',upper(mfilename)); end; if ~isnumeric(J) || ~isscalar(J) error('%s: "J" must be a scalar.',upper(mfilename)); end; if(J<1 || rem(J,1)~=0) error('%s: J must be a positive integer.',upper(mfilename)); end % Initialize the wavelet filters structure h = fwtinit(h,'ana'); %% ----- step 0 : Check inputs ------- definput.import = {'fwt'}; [flags,kv]=ltfatarghelper({},definput,varargin); %% ----- step 1 : Verify f and determine its length ------- % Change f to correct shape. % TO DO: if elements of h.a are not equal and the flag is 'per', % do zero padding to the next multiple of 2^J [f,Ls,W,wasrow,remembershape]=comp_sigreshape_pre(f,upper(mfilename),0); if(Ls<2) error('%s: Input signal seems not to be a vector of length > 1.',upper(mfilename)); end %% ----- step 2 : Check whether the input signal is long enough % input signal length is not restricted for expansive wavelet transform (extension type other than the default 'per') flen = length(h.filts{1}.h); if(strcmp(flags.ext,'per')) minLs = (h.a(1)^J-1)*(flen-1); % length of the longest equivalent filter -1 if Ls