exa130202.m

胡广书的《现代信号处理教程》附书光盘源代码。matlab实现。

% -----------------------------------------------------------------------－－
%   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