cwtnew.m

用于进行小波时频分析的matlab源码,可进行测井小波时频划分!

function coefs = cwtnew(signal,scales,wname,plotmode,xlim)
%CWT Real or Complex Continuous 1-D wavelet coefficients.
%   COEFS = CWT(S,SCALES,'wname') computes the continuous
%   wavelet coefficients of the vector S at real, positive
%   SCALES, using wavelet whose name is 'wname'.
%   The signal S is real, the wavelet can be real or complex. 
%
%   COEFS = CWT(S,SCALES,'wname','plot') computes
%   and, in addition, plots the continuous wavelet
%   transform coefficients.
%
%   COEFS = CWT(S,SCALES,'wname',PLOTMODE) computes and,
%   plots the continuous wavelet transform coefficients.
%   Coefficients are colored using PLOTMODE.
%   PLOTMODE = 'lvl' (By scale) or 
%   PLOTMODE = 'glb' (All scales) or
%   PLOTMODE = 'abslvl' or 'lvlabs' (Absolute value and By scale) or
%   PLOTMODE = 'absglb' or 'glbabs' (Absolute value and All scales)
%
%   CWT(...,'plot') is equivalent to CWT(...,'absglb')
%
%   You get 3-D plots (surfaces) using the same keywords listed
%   above for the PLOTMODE parameter, preceded by '3D'. For example:
%   COEFS = CWT(...,'3Dplot') or COEFS = CWT(...,'3Dlvl').
%
%   COEFS = CWT(S,SCALES,'wname',PLOTMODE,XLIM) computes, and
%   plots, the continuous wavelet transform coefficients.
%   Coefficients are colored using PLOTMODE and XLIM.
%   XLIM = [x1 x2] with 1 <= x1 < x2 <= length(S).
%
%   For each given scale a within the vector SCALES, the  
%   wavelet coefficients C(a,b) are computed for b = 1 to
%   ls = length(S), and are stored in COEFS(i,:)
%   if a = SCALES(i). 
%   Output argument COEFS is a la-by-ls matrix where la 
%   is the length of SCALES. COEFS is a real or complex matrix
%   depending on the wavelet type.
%
%   Examples of valid uses are:
%   t = linspace(-1,1,512);
%   s = 1-abs(t);
%   c = cwt(s,1:32,'cgau4');
%   c = cwt(s,[64 32 16:-2:2],'morl');
%   c = cwt(s,[3 18 12.9 7 1.5],'db2');
%   c = cwt(s,1:32,'sym2','lvl');
%   c = cwt(s,1:64,'sym4','abslvl',[100 400]);
%
%   See also WAVEDEC, WAVEFUN, WAVEINFO, WCODEMAT.

%   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.
%   Last Revision: 02-Feb-2000.
%   Copyright 1995-2002 The MathWorks, Inc.
%   $Revision: 1.18 $ $Date: 2002/03/28 17:24:50 $

% Check arguments.
if errargn(mfilename,nargin,[3:5],nargout,[0 1]), error('*'), end
err = 0;
if isempty(scales) ,         err = 1;
elseif min(size(scales))>1 , err = 1;
elseif min(scales)<eps,      err = 1;
end
if err
    errargt(mfilename,'Invalid Value for Scales !','msg');
    error('*')
end

precis = 10;
signal = signal(:)';
len    = length(signal);
coefs  = zeros(length(scales),len);
nbscales = length(scales);

[psi_integ,xval] = intwave(wname,precis);
wtype = wavemngr('type',wname);
if wtype==5 , psi_integ = conj(psi_integ); end

xval = xval-xval(1);
dx   = xval(2);
xmax = xval(end);
ind  = 1;
for k = 1:nbscales
    a = scales(k);
    j = [1+floor([0:a*xmax]/(a*dx))];
    if length(j)==1 , j = [1 1]; end
    f            = fliplr(psi_integ(j));
    coefs(ind,:) = -sqrt(a)*wkeep(diff(conv(signal,f)),len);
    ind          = ind+1;
end
if nargin<4 , return; end

dummyCoefs = coefs;
NBC = 128;
if strmatch('3D',plotmode)
    dim_plot = '3D';
else
    dim_plot = '2D';
end

if isequal(wtype,5)
   if ~isempty(findstr(plotmode,'lvl')) 
       plotmode = 'lvl';
   else
       plotmode = 'glb';   
   end
end
switch plotmode
  case {'lvl','3Dlvl'}
    lev_mode  = 'row';
    abs_mode  = 0;
    beg_title = ['By scale'];

  case {'glb','3Dglb'}
    lev_mode  = 'mat';
    abs_mode  = 0;
    beg_title = '';

  case {'abslvl','lvlabs','3Dabslvl','3Dlvlabs'}
    lev_mode  = 'row';
    abs_mode  = 1;
    beg_title = ['Abs. and by scale'];

  case {'absglb','glbabs','plot','2D','3Dabsglb','3Dglbabs','3Dplot','3D'}
    lev_mode  = 'mat';
    abs_mode  = 1;
    beg_title = ['Absolute'];

  otherwise
    plotmode  = 'absglb';
    lev_mode  = 'mat';
    abs_mode  = 1;
    beg_title = ['Absolute'];
    dim_plot  = '2D';
end

if abs_mode , dummyCoefs = abs(dummyCoefs); end
if nargin==5
    if xlim(2)<xlim(1) , xlim = xlim([2 1]); end    
    if xlim(1)<1   , xlim(1) = 1;   end
    if xlim(2)>len , xlim(2) = len; end
    indices = [xlim(1):xlim(2)];
    switch plotmode
      case {'glb','absglb'}
        cmin = min(min(dummyCoefs(:,indices)));
        cmax = max(max(dummyCoefs(:,indices)));
        dummyCoefs(dummyCoefs<cmin) = cmin;
        dummyCoefs(dummyCoefs>cmax) = cmax;

      case {'lvl','abslvl'}
        cmin = min((dummyCoefs(:,indices))')';
        cmax = max((dummyCoefs(:,indices))')';
        for k=1:nbscales
            ind = dummyCoefs(k,:)<cmin(k);
            dummyCoefs(k,ind) = cmin(k);
            ind = dummyCoefs(k,:)>cmax(k);
            dummyCoefs(k,ind) = cmax(k);
        end
    end
end

nb    = min(5,nbscales);
level = '';
for k=1:nb , level = [level ' '  num2str(scales(k))]; end
if nb<nbscales , level = [level ' ...']; end
nb     = ceil(nbscales/20);
tics   = 1:nb:nbscales;
tmp    = scales(1:nb:nb*length(tics));
labs   = num2str(tmp(:));
plotPARAMS = {NBC,lev_mode,abs_mode,tics,labs,''};

switch dim_plot
  case '2D'
    if wtype<5
	    titleSTR = [beg_title ' Coef a = ' level]; 
        plotPARAMS{6} = titleSTR;
        axeAct = gca;
        plotCOEFS(axeAct,dummyCoefs,plotPARAMS);
    else
        axeAct = subplot(2,2,1);
        titleSTR = sprintf('Real part of Ca,b for a = %s', level);
        plotPARAMS{6} = titleSTR;
        plotCOEFS(axeAct,real(dummyCoefs),plotPARAMS);
        axeAct = subplot(2,2,2);
        titleSTR = sprintf('Imaginary part of Ca,b for a = %s', level);
        plotPARAMS{6} = titleSTR;
        plotCOEFS(axeAct,imag(dummyCoefs),plotPARAMS);
        axeAct = subplot(2,2,3);
        titleSTR = sprintf('Modulus of Ca,b for a = %s', level);
        plotPARAMS{6} = titleSTR;
        plotCOEFS(axeAct,abs(dummyCoefs),plotPARAMS);
        axeAct = subplot(2,2,4);
        titleSTR = sprintf('Angle of Ca,b for a = %s', level);
        plotPARAMS{6} = titleSTR;
        plotCOEFS(axeAct,angle(dummyCoefs),plotPARAMS);
    end
    colormap(pink(NBC));

  case '3D'
    if wtype<5
	    titleSTR = [beg_title ' Values of Ca,b Coefficients for a = ' level]; 
        plotPARAMS{6} = titleSTR;
        axeAct = gca;
        surfCOEFS(axeAct,dummyCoefs,plotPARAMS);
        % xl = [1 len];
        % yl = [1 nbscales];
        % zl = [min(min(dummyCoefs)) max(max(dummyCoefs))];
        % set(axeAct,'Xlim',xl,'Ylim',yl,'Zlim',zl,'view',[-30 40]);
    else
        axeAct = subplot(2,2,1);
        titleSTR = sprintf('Real part of Ca,b for a = %s', level);
        plotPARAMS{6} = titleSTR;
        surfCOEFS(axeAct,real(dummyCoefs),plotPARAMS);
        axeAct = subplot(2,2,2);
        titleSTR = sprintf('Imaginary part of Ca,b for a = %s', level);
        plotPARAMS{6} = titleSTR;
        surfCOEFS(axeAct,imag(dummyCoefs),plotPARAMS);
        axeAct = subplot(2,2,3);
        titleSTR = sprintf('Modulus of Ca,b for a = %s', level);
        plotPARAMS{6} = titleSTR;
        surfCOEFS(axeAct,abs(dummyCoefs),plotPARAMS);
        axeAct = subplot(2,2,4);
        titleSTR = sprintf('Angle of Ca,b for a = %s', level);
        plotPARAMS{6} = titleSTR;
        surfCOEFS(axeAct,angle(dummyCoefs),plotPARAMS);
    end
end

%----------------------------------------------------------------------
function plotCOEFS(axeAct,coefs,plotPARAMS)

[NBC,lev_mode,abs_mode,tics,labs,titleSTR] = deal(plotPARAMS{:});

coefs = wcodemat(coefs,NBC,lev_mode,abs_mode);
img   = image(coefs');
%set(axeAct, ...
 %       'XTick',tics(1:4:length(tics)), ...
  %      'XTickLabel',labs(1:4:length(labs)), ...
   %     'XDir','normal', ...
    %    'Box','On' ...
     %   );
title(titleSTR,'Parent',axeAct);
ylabel('time b','Parent',axeAct);
xlabel('scales a','Parent',axeAct);
%----------------------------------------------------------------------
function surfCOEFS(axeAct,coefs,plotPARAMS)

[NBC,lev_mode,abs_mode,tics,labs,titleSTR] = deal(plotPARAMS{:});

img = surf(coefs);
set(axeAct, ...
        'YTick',tics, ...
        'YTickLabel',labs, ...
        'YDir','normal', ...
        'Box','On' ...
        );
title(titleSTR,'Parent',axeAct);
xlabel('time (or space) b','Parent',axeAct);
ylabel('scales a','Parent',axeAct);
zlabel('COEFS','Parent',axeAct);

xl = [1 size(coefs,2)];
yl = [1 size(coefs,1)];
zl = [min(min(coefs)) max(max(coefs))];
set(axeAct,'Xlim',xl,'Ylim',yl,'Zlim',zl,'view',[-30 40]);

colormap(pink(NBC));
shading('interp')
%----------------------------------------------------------------------