Wavelet analysis example from Matlab wavelet toolbox

Make plots for wavelet example from Matlab toolbox using electrical consumption time series leleccum (sampling rate 1/minute; tabulated over 2.7 days). This version does not require the wavelet toolbox, instead relying on Bretherton's own dwtHaar and idwtHaar functions.

Contents

Load data

  load leleccum;
  s = leleccum(1:4096);
  ls = length(s);

Plot signal

  figure(1)
  clf
  subplot(3,1,1)
  plot(s,'k')
  v = axis;
  v(1:2) = [0 ls];
  axis(v);
  ylabel('s')

Plot level-1 Haar transform of s

  % Single-level Haar (db1) wavelet transform

  ls = length(s);
  jo = 1:2:(ls-1);
  je = 2:2:ls;
  so = s(jo);
  se = s(je);
  cA1 = (so+se)/sqrt(2);
  cD1 = (so-se)/sqrt(2);
  C = [cA1 cD1];

  % Average and detail components of inverse transform
  A1 = zeros(1,ls);
  D1 = zeros(1,ls);
  A1(jo) = cA1/sqrt(2);
  A1(je) = cA1/sqrt(2);
  D1(jo) = cD1/sqrt(2);
  D1(je) = -cD1/sqrt(2);

  subplot(3,1,2)
  plot(1:ls/2,cA1,'b',(ls/2+1):ls,10*cD1,'r')
  v = axis;
  v(1:2) = [0 ls];
  axis(v);
  text(ls/4,0,'a^1')
  text(3*ls/4,500,'10d^1')
  ylabel('1-level Haar DWT')

Plot partitioning of signal into average and detail components.

  subplot(3,1,3)
  plot(1:ls,A1,'b',1:ls,10*D1,'r')
  v = axis;
  v(1:2) = [0 ls];
  axis(v);
   text(ls/2,750,'A^1')
  text(ls/2,-250,'10D^1')
  ylabel('Partitioning of signal')

3-level Haar transform

  maxlev = 3;
  C = dwtHaar(s,maxlev);

  D = idwtHaar(C,maxlev);

  % Extract information at each scale. Signal s = sum(D)

  D1 = D(1,:);
  D2 = D(2,:);
  D3 = D(3,:);
  A3 = D(4,:);

Plot signal and the 3rd level averaged signal

  figure(2)
  clf
  subplot(3,1,1)
  plot(1:ls,s,'k',1:ls,A3,'b')
  v = axis;
  v(1:2) = [0 ls];
  axis(v);
  legend('signal','A^3')

  % Plot level-3 Haar transform

  subplot(3,1,2)
  plot(1:ls/8,C(1:ls/8),'b',(ls/8+1):(ls/4),10*C((ls/8+1):(ls/4)),'c',...
       (ls/4+1):(ls/2),10*C((ls/4+1):(ls/2)),'m',...
       (ls/2+1):ls,10*C((ls/2+1):ls),'r')
  v = axis;
  v(1:2) = [0 ls];
  axis(v);
   ylabel('3-level Haar DWT')
  legend('a^3','10d^3','10d^2','10d^1')

  % Plot the three detail components of signal

  subplot(3,1,3)
  yoff = 20;
  plot(1:ls,D1,'r',1:ls,D2+yoff,'m',1:ls,D3+2*yoff,'c')
  v = axis;
  v(1:2) = [0 ls];
  axis(v);
  axis([0,ls,-20,60])
  legend('D1','D2','D3')

Signal compression

  % Calculate complete p-level Haar transform

  P = log2(ls);
  C = dwtHaar(s,P);

  % Plot the largest nkeep wavelet coefficients, for various
  % choices nkeep

  figure(3)
  clf
  subplot(2,2,1)
  Csq = C.^2;
  [mCsqsort,Isort] = sort(-Csq);  % Tricky way to sort wavelet coeffs from
                                % largest to smallest in magnitude.
  cumpow = cumsum(mCsqsort)/sum(mCsqsort); % Cum frac of Wavelet power
  nkeep = [10 30 100];
  semilogy(1:nkeep(1),1-cumpow(1:nkeep(1)),'r',...
           nkeep(1):nkeep(2),1-cumpow(nkeep(1):nkeep(2)),'g',...
           nkeep(2):nkeep(3),1-cumpow(nkeep(2):nkeep(3)),'b');
  xlabel('Retained wavelet coeffs.')
  ylabel('Resid. Power Frac.')
  legend(int2str(nkeep(1)),int2str(nkeep(2)),int2str(nkeep(3)))

Reconstruct signal with all but largest nkeep wavelet coeffs zeroed out

and plot the compression error

  subplot(2,2,2)
  col = 'rgb';  % Vector of plotting colors
  for ik = 1:length(nkeep)
    jsort = Isort(1:nkeep(ik));
    C2 = C - C;  % Zero wavelet coeffs of compressed signal
    C2(jsort) = C(jsort); % ...except for the nkeep(ik) largest.
    s2 = sum(idwtHaar(C2,P)); % Reconstruction after compression

    % Plot compression error
    plot(1:ls, s2-s, col(ik))
    if (ik==1)
      v = axis;
      v(1:2) = [0 ls];
      axis(v);
      xlabel('Time index [min]')
      ylabel('Compression error')
      hold on
    end
  end
%  legend(int2str(nkeep(1)),int2str(nkeep(2)),int2str(nkeep(3)))
  hold off

  % Plot levels of largest wavelet coeffs
  subplot(2,1,2)
  loglog(1:ls,Csq,'k',Isort(1:nkeep(3)),-mCsqsort(1:nkeep(3)),'bx',...
           Isort(1:nkeep(2)),-mCsqsort(1:nkeep(2)),'gx',...
           Isort(1:nkeep(1)),-mCsqsort(1:nkeep(1)),'rx')
  v = axis;
  v(1:2) = [0 ls];
  axis(v);
  ylabel('squared level-12 wavelet coeff.')
  xlabel('Index in wavelet coeff vector C')
  legend('No trunc',int2str(nkeep(3)),int2str(nkeep(2)),int2str(nkeep(1)))
  hold on
  ytxt = v(3)^0.9 * v(4)^0.1;
  text(1.1,ytxt,'a^{12}')
  for lev = 1:P
     plot([2^(P-lev)+0.5 2^(P-lev)+0.5],[v(3) v(4)],'k--')
     text(1.55*2^(P-lev),ytxt,['d^{' num2str(lev) '}'])
  end
  hold off

Published with MATLAB® R2012b