% -----------------------------------------------------------------------－－
%   exa130202.m,   例13.2.2
%  求沖激函數、階躍函數及三角函數的小波變換、模極大線并求它們的李氏指數
%  注：在該程序中，用到了子程序 RWT.m，
%      該程序請讀者在如下的網站上下載:
%                      http://www-stat.stanford.edu/~wavelab/
%      因為該網站上的小波分析軟件并沒有列入MATLAB的工具箱，涉及到知識產權問題，
%      因此不能將其列入國內公開出版物上，但讀者可以自由下載。
%------------------------------------------------------------------------
close all;

nvoice = 16;
%   nvoice表示每2倍s之間被分為多少點

wavelet = 'Sombrero'; 
% 選擇小波：可選擇的小波包括： 'Gauss','DerGauss','Morlet'and'Sombrero'(mexh)
    N = 64*8;
    n = N;
    oct = 5;
    scale = 8;
   	noctave = floor(log2(n))-oct;
	nscale  = nvoice .* noctave;
    ytix   = linspace(2+(oct-floor(log2(scale))),log2(N)+2-floor(log2(scale)),nscale);
%  ytix 從imageRWT函數得到 理論上應等于log2(s)
%  ytix = linspace(log2(scale)-noctave,log2(scale),nscale);
    xtix   = linspace(0,N,N);
for count=1:3
   if count==1
        sig = zeros(1,N);sig(N/2) =sig(N/2) +1;% 沖激函數
        ti='沖激函數';
    elseif count==2
        sig = zeros(1,N);sig(N/4:3*N/4) = sig(N/4:3*N/4)+1; % 矩形窗
        ti='矩形窗';
    else
        sig(N/4:N/2) = linspace(0,1,N/4+1);
        sig(N/2:3*N/4) = linspace(1,0,N/4+1); % 三角窗
        ti='三角窗';
    end
    sig = sig+1e-8*randn(1,N);
    rwt = RWT(sig,nvoice,wavelet,oct,scale);
    [n,nscale] = size(rwt);
   
    subplot(4,3,count); % 顯示信號波形
    plot(sig);grid;title(ti);
    axis([1 N min(sig)-.1 max(sig)+.1]);
    subplot(4,3,3+count); % 顯示小波變換結果
    imagesc(xtix,ytix,flipud(rwt'));colormap(gray);
     maxmap = MM_RWT(rwt);
%  maxmap就是與rwt相同大小的矩陣，但在maxima處為1，其它處為0

   [skellist,skelptr,skellen] = SkelMap(maxmap);
% skellist 依次記錄了每條"山脊"點坐標對(j,i):[j1,i1,j2,i2,...]。
% 則所有最大點的坐標都存在里面
% skelptr 為指向skellist的"指針"skelptr第n個元素是指第n條"山脊"的起點在skellist中的位置。
% 如skelptr(2)=65，skelptr(3)=129，則skellist(65:128)就是第2條"山脊"上各點的坐標對。
% skellen(n)就是第n條"山脊"的在skellist中的長度。
 
    block=rwt';
%   找出最大值及其位置
    [maxvalue,maxpos]=max(block,[],2);
%   找出最小值及其位置
    [minvalue,minpos]=min(block,[],2);
    subplot(4,3,6+count); % 繪制最大最小線
    PlotSkelMap(n,nscale,skellist,skelptr(:),skellen(:),'','b',[],nvoice,min(ytix),noctave);grid;
    ylabel('')
%    ylabel('log2(s)');
    subplot(4,3,9+count);% 估計斜率
%    plot(ytix,log2(abs(flipud(minvalue)))-ytix'/2,'b-');hold on;
%   藍色實線為最小線，  -ytix'/2這項相當于除以sqrt(s)
%   
    plot(ytix,log2(abs(flipud(minvalue))),'b-');hold on;

%      在log2()之內./sqrt(2.^ytix')與在log2()之外-ytix'/2等效
    plot(ytix,log2(abs(flipud(maxvalue))),'r--');
    grid on;
    hold off;
 end