scalwav.m

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

function [s,w] = scalwav(lpf,iter)
%Orthonormal scaling function and wavelet generating routine.
%
%[s,w] = scalwav(lpf,iter) takes the lowpass filter (lpf)
%of an orthonormal filter bank, and returns samples of the scaling 
%function and the corresponding wavelet. The scaling 
%function and wavelet samples are returned in arrays 
%s and w respectively.
%
%The variables used are:
%lpf : the low pass filter. It must be a row vector. 
%iter: number of iterations to be performed. The greater 
%      number of iterations will give a greater number of samples
%      of scaling function and wavelet.
%
%Example: 
%        Suppose lpf =[0.4830, 0.8365, 0.2241, -0.1294]
%
%Then typing 
%[s,w] scalwav(lpf,10)
%
%at the MATLAB prompt will compute the first 10 iterations of the 
%D-4 scaling function and wavelet and return the computed samples of 
%the scaling function and wavelet in the arrays, s and w, respectively. 
%
%Refer to Chapter 3 for information on orthonormal scaling 
%functions and wavelets. 
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%

   e = lpf*lpf';
   lpf = lpf/sqrt(e); %normalize the energy for convergence
   l = length(lpf)-1; %length of the scaling 
                      %function and wavelet in 
                      %time

    for i=0:l
      hpf(i+1) = (-1)^i*lpf(l+1-i);
    end; %endfor ...giving the highpass filter

   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 of the 
                              %wavelet function.

   %remaining iterations performed here
   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; %construction is complete

  len = length(s);
  s =sqrt(2)*[s 0]; 
  w =sqrt(2)*[w 0];
  x = 0:l/len:l; %the x coordinates of the 
                       %scaling function and the 
                       %wavelet samples computed

%plot the scaling function  
  figure(1);plot(x,s);title('Scaling Function');xlabel('t');
%plot the wavelet function 
  figure(2);plot(x,w);title('Wavelet');xlabel('t');