Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Diff of /wavelets/fwt.m [000000] .. [77a22a] Maximize Restore

  Switch to side-by-side view

--- a
+++ b/wavelets/fwt.m
@@ -0,0 +1,157 @@
+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<minLs
+     error('%s: Input signal length is %d. Minimum signal length is %d or use %s flag instead. \n',upper(mfilename),Ls,minLs,'''ppd''');
+   end;
+end
+
+
+%% ----- step 3 : Run computation
+ c = comp_fwt_all(f,h.filts,J,h.a,'dec',flags.ext);
+
+