biscwav.m

这是伯克里wavelet transforms一书中的例子的代码

function [s,w] = biscwav(lpf,hpf,iter)
%Biorthogonal scaling function and wavelet generating routine.
%
%[s,w] = biscwav(lpf,hpf,iter) takes the low pass and the high 
%pass filters of a biorthogonal filter bank and returns samples of 
%the corresponding scaling function and wavelet in arrays s and w
%respectively.
%
%The variables used are
%lpf : Low pass filter of the biorthogonal filter bank.
%hpf : High pass filter of the biorthogonal filter bank.
%iter: Number of iterations to be performed. Greater number of
%      iterations will give a greater number of points for the 
%      the scaling function and wavelet.
%
%Example:
%Suppose
%lpf = [0.0378, -0.0238, -0.1106, 0.3774, 0.8527, 0.3774, -0.1106,
%       -0.0238, 0.0378] and
%hpf = [0.0645, -0.0407, -0.4181,  0.7885, -0.4181, -0.0407, 0.0645]
%
%Then, 
%[s,w] = biscwav(lpf,hpf,10)
%
%computes the first 10 iterations of the scaling function and wavelet
%corresponding to the given lpf and hpf.
%
%Refer to Chapter 4 for information on biorthogonal scaling 
%functions and wavelets.
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%

  ls = length(lpf)-1; %support of the scaling function
  lw = (ls +length(hpf)-1)/2; %support of the wavelet function

    el = lpf*lpf';
    eh = hpf*hpf';
    lpf = lpf/sqrt(el);
    hpf = hpf/sqrt(eh);

    l1 = 2*length(lpf)-1;
%upsample
    s1 = zeros([1 l1]);
    s1(1:2:l1) = lpf;
  s = sqrt(2)*conv(s1,lpf); %first iteration of the scaling function

   l1 = 2*length(hpf)-1;
%upsample
   w1 = zeros([1 l1]);
   w1(1:2:l1) = hpf;
   w = sqrt(2)*conv(w1,lpf); %first iteration for the wavelet

 %begin rest of the iterations
  for i=1:(iter-1)
    l1 = 2*length(s)-1;
%upsample
    s1 = zeros([1 l1]);
    s1(1:2:l1) = s;
    s  = sqrt(2)*conv(s1,lpf); %scaling function

    l1 = 2*length(w)-1;
%upsample
    w1 = zeros([1 l1]);
    w1(1:2:l1) = w;
    w = sqrt(2)*conv(w1,lpf); %wavelet function
 
  end %end for

  l1 = length(s); 
  l2 = length(w);
  s = sqrt(el)*[s 0];
  w = sqrt(eh)*[w 0];

  xs = 0:ls/l1:ls; %support for scaling function
  xw = 0:lw/l2:lw; %support for wavelet function
 %Plot the scaling function
 figure(1);plot(xs,s);xlabel('t');title('Scaling function');
%...and the wavelet
 figure(2);plot(xw,w);xlabel('t');title('Wavelet');