function [mi,wi] = meyer(len,epsi)
%Meyer scaling function and wavelet generating routine.
%
%[m,w] = meyer(length,epsi) produces two arrays, m and w 
%that contain samples of a Meyer type scaling function and
%wavelet respectively. The shape of the scaling function and
%the wavelet are defined by the variable epsi. The power spectrum
%of the scaling function is of the raised cosine type. 
%
%The equation of the transition band of the power spectrum of the 
%scaling function is:
%
%|PHI(w)|^2 =  (1 + cos((w - pi + 2*pi*epsi)/(4*epsi)))/2
%pi - 2*pi*epsi <= w <= pi + 2*pi*epsi
%
%The input variables are:
%length : The desired length of the scaling function and wavelet.
%         The length MUST be even.
%epsi   : A variable that can take values between 0 and 1/6. It
%         determines the width of the transition band of the 
%         Meyer type scaling function power spectrum. 
%
%Example:
%[m,w] = meyer(256,1/12)
%
%computes 256 samples of a Meyer type scaling function and wavelet. The 
%transition band of the scaling function is from pi-2*pi*(1/12)=5*pi/6
%to pi+2*pi*(1/12)=7*pi/6.
%
%Refer to Chapter 3 for information on the Meyer type scaling
%function and wavelet.
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%

  if (len/2 ~= round(len/2)) %if length is odd then
     len = len +1; %round it to be even-one greater than the input odd
                   %length.
     fprintf('Odd length input: using length = %3d instead\n',len);
  end %endif

  if (epsi > 1/6)
     epsi = 1/6;
     fprintf('epsi value greater than 1/6, using default value of 1/6 instead\n');
  elseif (epsi < 0)
    epsi = 0;
    fprintf('epsi value less than 0, using default value of 0 instead\n');
  end %endif

  domega = 4*pi/len; %fix the frequency increments

  for i = 1:len/2 %you need to compute only half of the filter
   
   incr = (i-1)*domega;
   %stopband
   if (incr > pi +2*pi*epsi)
     m(i) = 0;

   %Raised cosine type transition band
   elseif (incr>pi-2*pi*epsi & incr <= pi +2*pi*epsi)
     m(i) = (1+cos((incr - pi+2*pi*epsi)/(4*epsi)))/2;

   %if the eps variable is specified to be zero
   else
     if ((incr == pi) & (epsi == 0))
       m(i) = 0.5;
     else   
       m(i) =1;
     end %end inner if

   end %endif
  end %endfor

   m = sqrt([m 0 m(len/2:-1:2)]); %this gives the Meyer Scaling 
                                     %function spectrum

   mr =rright(m,len/2); %the magnitude spectrum of the HPF
   
   ws = mr(1:2:len).*exp(-sqrt(-1)*domega*(0:2:len-1)/2);
   %HPF Fourier Transform

   w = [ws ws] .* m; %wavelet spectrum

   mi = real(ifft(m)); %scaling function
   wi = real(ifft(w)); %wavelet function

   mi = rright(mi,len/2);%rotate right to center 
   wi = rright(wi,len/2 -2);%the two sequences

%for smooth plots

   if(len<=8) %for shorter lengths...
     m = interp(mi,16,1,2/3);
     w = interp(wi,16,1,2/3);
   elseif(len > 8 & len <= 64) %for longer lengths
     m = interp(mi,4);
     w = interp(wi,4);
   else
     m = mi;
     w = wi;
   end %endif...
   l = length(m); %length for plotting..
   x = -len/2:len/(l):(len/2 - len/(2*l));
  
   figure(1);plot(x,m);title('Meyer Scaling Function'); %Plot the scaling function and wavelet
   figure(2);plot(x,w);title('Meyer Wavelet');