Scilab Wavelet Toolbox

Scilab Wavelet Toolbox Git Source Tree

Root/macros/wdencmp.sci

1function [xc,cxc,lxc,perf0,perfl2] = wdencmp(o,varargin)
2// De-noising or compression using wavelets.
3// Calling Sequence
4// [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('gbl',X,'wname',N,THR,SORH,KEEPAPP)
5// [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('gbl',C,L,W,N,THR,SORH,KEEPAPP)
6// [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('lvd',X, 'wname',N,THR,SORH)
7// [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('lvd',C,L, 'wname',N,THR,SORH)
8// [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('lvd',X, 'wname',N,THR,SORH)
9// [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('lvd',C,L, 'wname',N,THR,SORH)
10// Parameters
11// X : input signal (1-D or 2-D)
12// C,L:wavelet decomposition structure
13// THR : positive threshold
14// XC: de-noised or compressed version of X
15// [CXC,LXC]: wavelet decomposition of XC
16// PERFL2 and PERF0: are L^2 recovery and compression scores in percentages.
17// PERFL2: PERFL2 = 100*(vector-norm of CXC/vector-norm of C)^2
18// N: level of Wavelet decomposition
19// 'wname': is a string containing the wavelet name.
20// SORH('s' or 'h') : soft or hard thresholding
21// KEEPAPP: = 1, approximation coefficients cannot be thresholded, otherwise it is possible.
22// 'gbl': using one threshold value
23// 'lvd': level-dependent thresholds (THR must be of length N). For 2-D case THR must be a matrix of size 3 by N in the three orientations horizontal, diagonal and vertical.
24//
25// Description
26// performs a de-noising or compression process
27// of a signal or an image using wavelets.
28//
29// Examples
30// x=sin(2*%pi*(0:0.01:1));
31// xn=x+rand(x);
32// thr=35;
33// //compression using global thresholding
34// [xcomp,cxd,lxd,perf0,perfl2] = wdencmp('gbl',xn,'db3',3,thr,'h',1);
35// //denoising
36// // Find default values
37// [thr,sorh,keepapp] = ddencmp('den','wv',x);
38// // De-noise signal using global thresholding option.
39// xd = wdencmp('gbl',xn,'db3',4,thr,sorh,keepapp);
40// scf();clf();
41// plot(xn);
42// plot(xcomp,'r');
43// plot(xd,'g');
44// legend(["noisy signal","compressed signal","de-noised signal"],1);
45// See also
46// ddencmp
47// wavedec
48// wavedec2
49// wden
50// wthresh
51// Authors
52// H. Nahrstaedt - 2010-2012
53
54
55
56[nargout,nargin]=argn();
57dim = 1; // initialize dimension to 1D.
58nbIn = nargin;
59nbOut = nargout;
60select o
61 case 'gbl' , minIn = 7; maxIn = 8;
62 case 'lvd' , minIn = 6; maxIn = 7;
63 else
64 error('Invalid argument value.')
65end
66
67okOut = [0:1 3:5];
68if ~or(okOut==nbOut)
69 error('Invalid number of output arguments.');
70end
71
72if nbIn == minIn
73 x = varargin(1); indarg = 2;
74 if min(size(x))~=1, dim = 2; end
75else
76 c = varargin(1); l = varargin(2); indarg = 3;
77 if min(size(l))~=1, dim = 2; end
78end
79
80// Get Inputs
81w = varargin(indarg);
82n = varargin(indarg+1);
83thr = varargin(indarg+2);
84sorh = varargin(indarg+3);
85if (o=='gbl') , keepapp = varargin(indarg+4); end
86
87
88// Wavelet decomposition of x (if not given).
89if ((o=='gbl') & nbIn==7) | ((o=='lvd') & nbIn==6)
90 if dim == 1, [c,l] = wavedec(x,n,w);
91 else [c,l] = wavedec2(x,n,w);
92 end
93end
94
95// Wavelet coefficients thresholding.
96if (o=='gbl')
97 if keepapp
98 // keep approximation.
99 cxc = c;
100 if dim == 1, inddet = l(1)+1:length(c);
101 else inddet = prod(l(1,:))+1:length(c); end
102 // threshold detail coefficients.
103 cxc(inddet) = wthresh(c(inddet),sorh,thr);
104 else
105 // threshold all coefficients.
106 cxc = wthresh(c,sorh,thr);
107 end
108else
109 if dim == 1, cxc = wthcoef('t',c,l,1:n,thr,sorh);
110 else
111 cxc = wthcoef2('h',c,l,1:n,thr(1,:),sorh);
112 cxc = wthcoef2('d',cxc,l,1:n,thr(2,:),sorh);
113 cxc = wthcoef2('v',cxc,l,1:n,thr(3,:),sorh);
114 end
115end
116lxc = l;
117
118// Wavelet reconstruction of xd.
119if dim == 1,xc = waverec(cxc,lxc,w);
120else xc = waverec2(cxc,lxc,w);
121end
122
123if nbOut<4 , return; end
124
125// Compute compression score.
126perf0 = 100*(length(find(cxc==0))/length(cxc));
127if nbOut<5 , return; end
128
129// Compute L^2 recovery score.
130nc = norm(c);
131if nc<%eps
132 perfl2 = 100;
133else
134 perfl2 = 100*((norm(cxc)/nc)^2);
135end
136endfunction

Archive Download this file

Branches