smallexp.m

在matlab的小波分析中标,设计一个包含若干不连续点的一维信号

function [ Lips ] = smallExp( )
% =========================================
%                  习题8.2
%       计算信号奇异点的Lipschitz指数
% =========================================

close all;

% 构造信号
% 引用wavelab中的数据,即A Wavelet Tour of Signal Processing(2nd edition)中fig6.6的信号
sig = GetSignal;
n = 512;        % 信号长度

% % 引用实际图像信号
% img = imread('lenna.jpg');
% sig = double(img(30,:));
% n = 256;        % 信号长度

% 计算信号在多个尺度下的小波变换wt
s_scale = 1;    % 最小尺度
l_scale = 32;   % 最大尺度
wt = CWT(sig,s_scale:l_scale,'gaus2');
   
% 计算模极大点
maxmap = MM_WT(wt,1);
% 连接模极大曲线
[skellist,skelptr,skellen] = SkelMap(maxmap);

% 图1：信号图和小波变换图
% 信号图和小波变换图的位置
delta = 1/15;
unit = (1-3*delta)/3;
h1 = [delta 2*(unit+delta) 1-2*delta unit];
h2 = [delta delta 1-2*delta 2*unit];

figure(1);
clf;
% 画出信号
axes('position',h1);    
plot(sig);
axis([1 n min(sig)-.1 max(sig)+.1]);
ylabel('f(t)');
% 画出信号的小波变换结果
axes('position',h2);
ImageWT(wt,n,s_scale,l_scale);

% 图2，画出模极大曲线
figure(2);
axis('ij');
xlabel('u');
ylabel('log2(s)');
plotsymb = ['k' '-'];
hold on;
nchain  = length(skelptr);
for k=1:nchain
	ix  = skelptr(k): (skelptr(k) + skellen(k)-1);
	vec = skellist(:,ix);
    plot(vec(2,:),log2(vec(1,:)),plotsymb);
%	plot(vec(2,:),log2(vec(1,:))-log2(n),plotsymb); % 注意纵坐标,转化回[0 1]之间的信号的变换尺度
end	
hold off;
    
% 图3，绘制极大曲线在log2(s),log2|WT(s,u)|平面上的曲线图  
figure(3);
xlabel('log2(s)');
ylabel('log2|WT(s,u)|');
hold on;
for k=1:nchain
    ix  = skelptr(k): (skelptr(k) + skellen(k)-1);
	vec = skellist(:,ix);
    pvec = abs( wt((vec(2,:)-1)*32 + vec(1,:)) );
    %plot(log2(vec(1,:))-log2(n),log2(pvec),plotsymb);
    if(k==5)
        plot(log2(vec(1,:)),log2(pvec),plotsymb);
        % plot(log2(vec(1,:))-log2(n),log2(pvec),plotsymb);   
    end
    if(k==6)
        plot(log2(vec(1,:)),log2(pvec),['k' '--']);  
        % plot(log2(vec(1,:))-log2(n),log2(pvec),['k' '--']);   
    end
    % 计算Lipschitz指数
    if(skellen(k)<=4)
        m1 = 1;
        m2 = min(4,skellen(k));
    else
        m1 = 4;
        m2 = min(16,skellen(k));
    end
    [P S] = polyfit(log2(vec(1,m1:m2)),log2(pvec(m1:m2)),1);
    % [P S] = polyfit(log2(vec(1,:)),log2(pvec),1);
    Lips(1:2,k) = [skellist(2,skelptr(k)) ; P(1)-0.5];
end   
% legend('5','6',-1);
hold off;