Lipschitz指数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;