deno_uv.m

Matlab toolbox wave 3example小波计算,希望对大家有用

% This function de-noises an input of length N = 2^n using 
% the wavelet shrinkage method in which the "universal" 
% threshold is used.
% Note: The function requires the use of UVi_Wave wavelet toolbox.
% Inputs: x0 -- Original signal.
%              r -- White noise.
%            h0 -- Impulse response of the lowpass analysis filter.
%            h1 -- Impulse response of the highpass analysis filter.
%             f0 -- Impulse response of the lowpass synthesis filter.
%             f1 -- Impulse response of the highpass synthesis filter.
%  choice:  0 -- Hyperbolic (almost hard) shrinkage.
%              1 -- Donoho's soft shrinkage.
%  Output: y -- Denoised signal. 
% Written by W.-S. Lu, University of Victoria.
% Last modified March 1, 2003.

function y = deno_uv(x0,r,h0,h1,f0,f1,choice)

x0 = x0(:);
r = r(:);
x = x0 + r;
Nx = length(x);
t = 0:1/(Nx-1):1;

% Handle the not-power-of-2 case

p = ceil(log10(Nx)/log10(2));
N = 2^p;
if N > Nx,
 xa = [x;zeros(N-Nx,1)];
 else xa = x;
end

% Wavelet decomposition

xd = wt(xa,h0,h1,p);

% Estimate noise variance

d1 = xd(N/2+1:N);
db = mean(d1);
s2 = sum((d1-db).^2)/(N/2-1);

% Compute Donoho's universal threshold

delta = sqrt(2*log10(N)*s2)

% Soft (or hyperbolic) shrinkage of wavelet coefficients xd

if choice == 1,
 xs = sign(xd).*max(abs(xd)-delta,zeros(N,1));
else
 xs = sign(xd).*sqrt(max(xd.^2-delta^2,zeros(N,1)));
end

% Reconstruct the denoised signal

y1 = iwt(xs,f0,f1,p);
y = y1(1:Nx);
y = y(:);

% Performance evaluation

nf = norm(x0);
disp('signal-to-noise ratio before denoising:')
SNR_before = 20*log10(nf/norm(r))
disp('signal-to-noise ratio after denoising:')
SNR_after = 20*log10(nf/norm(y-x0))

% Plot the results

mi = min([min(x0) min(x) min(y)]);
mi = (1 - 0.1*sign(mi))*mi;
ma = max([max(x0) max(x) max(y)]);
ma = (1 + 0.1*sign(ma))*ma;
figure(1)
subplot(221)
plot(t,x0)
title('original signal')
axis([0 1 mi ma])
axis('square')
grid
subplot(222)
plot(t,r)
title('the noise')
axis('square')
grid
subplot(223)
plot(t,x)
title('noise-contaminated signal')
axis([0 1 mi ma])
axis('square')
grid
subplot(224)
plot(t,y)
title('denoised signal')
axis([0 1 mi ma])
axis('square')
grid