cdeowavf(varargin);

function varargout = cdeowavf(varargin);
% cdeo  wavelet -- TYPE 5 WAVELET according to waveinfo
%   [PSI,T] = cdeo(LOWB,UPPB,N,alpha) returns deO
%   complex wavelet function evaluated on an
%   N point regular grid on the interval [LOWB,UPPB].
%   Output arguments are the wavelet function PSI
%   computed on the grid T (defined by LOWB, UPPB and N)
%
%   These functions have [-12 12] as effective support.
%
%   N must be a power of two.
%   See also CDEOINFO, WAVEFUN, WAVEINFO.
%
%   Luciana R. Soares, H.M. de Oliveira.
%   Last Revision: 10-Dec-2002

% Check arguments
if errargn(mfilename,nargin,[4],nargout,[0:2]), error('*'), end

% Translate input arguments...

    [lowb,uppb,N,label] = deal(varargin{:});
    ind = strncmpi('cdeo',label,4);
    if isequal(ind,1) label([1:4]) = []; end
    alpha = wstr2num(deblank(label));
 
% Check alpha factor (roll-off)
if isempty(alpha) error('** deO: Invalid wavelet number!'); end
if (alpha < 0) | (alpha > 1/3)
    error('** deO: Invalid value for alpha **')
end
if (alpha == 0)
   alpha = 10.^-5; % alpha ~ 0
   disp('deO: Shannon square root Wavelet');
end

% Check N
if errargt(mfilename,log(N)/log(2),'int')
    error('** deO: Invalid value for N **')
end

% Check Interval
if errargt(mfilename,uppb-lowb,'re0')
   error('** deO: Verify upper/lower bound... **')
end

% Transform interval bounds to grid
lint = (uppb-lowb);
step = lint/N;
if step > 0.3
   disp('** deO: the number of grid points is not enough for that support... Be carefull **');
end
t  = [lowb:step:uppb];
td = t - 1/2;

% Initial Variables
va = 2*(1 - alpha);
vb = 1 + alpha;
vc = 1 - alpha;
vd = 2*(1 + alpha);
epslon = 0.0001;

% Computing critical values...

crit_cb1 = 1/(2*(vc-vb));
crit_cb2 = -1/(2*(vc-vb));

y = find(td==0);
ycb1 = find(abs(td - crit_cb1) < epslon);
ycb2 = find(abs(td - crit_cb2) < epslon);

crit_da1 = 1/(2*(vd-va));
yda1 = find(abs(td - crit_da1) < epslon);
crit_da2 = -1/(2*(vd-va));
yda2 = find(abs(td - crit_da2) < epslon);

% Avoiding critical values...

td(y) = epslon;
td(ycb1) = crit_cb1 + epslon;
td(ycb2) = crit_cb2 + epslon;
td(yda1) = crit_da1 + epslon;
td(yda2) = crit_da2 + epslon;

% Wavelet function (psi)

   % Real part
   hra = sin(va*pi*td)./(pi*td);
   hrb = sin(vb*pi*td)./(pi*td);
   mrbc = (2*abs(vc-vb)*(cos(vc*pi*td) + 2*(vc-vb).*td.*sin(vb*pi*td)))./(pi.*(1-(2*td*(vc-vb)).^2));
   mrad = (2*abs(vd-va)*(cos(vd*pi*td) + 2*(vd-va).*td.*sin(va*pi*td)))./(pi.*(1-(2*td*(vd-va)).^2));
   RdeO = (hra - hrb + mrbc + mrad)/(2*sqrt(2*pi));
   % Imaginary part
   denom = pi * td;
   hia = cos(va*pi*td)./denom;
   hib = cos(vb*pi*td)./denom;
   den1 = pi.*(1-(2*td*(vc-vb)).^2);
   den2 = pi.*(1-(2*td*(vd-va)).^2);
   mibc = (2*abs(vc-vb)*(sin(vc*pi*td) - 2*(vc-vb).*td.*cos(vb*pi*td)))./(den1);
   miad = (2*abs(vd-va)*(sin(vd*pi*td) - 2*(vd-va).*td.*cos(va*pi*td)))./(den2);
   IdeO = (hia - hib + mibc + miad)/(2*sqrt(2*pi));
   % deO Wavelet
   psi = 2*sqrt(pi)*(RdeO + j*IdeO);
end

% Set output arguments
if nargout ~= 2 
   error('**deO: too few output parameters: should be [psi, t]=...'); 
end
varargout = {psi,t};
   

%% typical deo computation
%% load a signal (e.g. cuspamax, nelec, lelecum etc.)
% load cuspamax;
% dados = cuspamax(1:1024);
% lv = length(dados);
% subplot(311), plot(dados); title('signal...');
%% performing a cwt by deo0.33333
%% all integer scales from 1 to 64
% subplot(313)
% wt = cwt(dados,1:64,'deo0.33333','plot');
% title('CWT, absolute coefficients')
% colormap(bone(64));
% ylabel('scale');