smoothn

Fast and robust smoothing of one-dimensional and multidimensional data

Download
Syntax
Description
Important note
1-D smoothing
2-D smoothing
2-D smoothing with missing data
3-D smoothing
Cardioid
Reference
See also
About the author

Download

smoothn.m

Syntax

   z = smoothn(y)
   z = smoothn(y,s)
   z = smoothn(y,W)
   z = smoothn(y,W,s)
   z = smoothn(...,'robust')
   z = smoothn(...,'option_name1',option_value1,...)
   [z,s,exitflag] = smoothn(...)

Description

z = smoothn(y) automatically smoothes the uniformly smpled array y. y can be any N-D noisy array (time series, images, 3D data,...). Non finite data (NaN or Inf) are treated as missing values.

z = smoothn(x,s) smoothes the array y using the smoothness parameter s. Parameter s must be a real positive scalar. The larger s is, the smoother the output will be. If s is omitted (see previous syntax) or empty (i.e. s = []), it is automatically determined using the generalized cross-validation (GCV) method.

z = smoothn(y,W) or z = smoothn(y,W,s) specifies a weighting array W of real positive values, that must have the same size as y. Note that a nil weight corresponds to a missing value.

z = smoothn(...,'robust') carries out a robust smoothing that minimizes the influence of outlying data.

[z,s] = smoothn(...) also returns the calculated value for s so that you can fine-tune the smoothing subsequently if needed.

An iteration process is used in the presence of weighted and/or missing values. z = smoothn(...,Option_Name,Option_Value) smoothes with the termination parameters specified by Option_Name and Option_Value. They can contain two criteria:

TolZ: Termination tolerance on Z (default = 1e-3), TolZ must be in ]0,1[
MaxIter: Maximum number of iterations allowed (default = 30)

Syntax: [z,...] = smoothn(...,'MaxIter',50,'TolZ',1e-4)

[z,s,exitflag] = smoothn(...) returns a boolean value exitflag that describes the exit condition of smoothn:

exitflag = 1: smoothn converged.
exitflag = 0: Maximum number of iterations was reached.

Class Support: Input array can be numeric or logical. The returned array is of class double.

Important note

The N-D (inverse) discrete cosine transform functions (dctn.m and idctn.m) are required.

1-D smoothing

Create noisy 1-D data and smooth them using smoothn. The smoothness parameter s is determined automatically. The plot shown below compares regular smoothing with a robust smoothing using the 'robust' option.

x = linspace(0,100,2^8);
y = cos(x/10)+(x/50).^2 + randn(size(x))/5;
y([70 75 80]) = [5.5 5 6];
[z,s] = smoothn(y); % Regular smoothing
zr = smoothn(y,'robust'); % Robust smoothing
subplot(121), plot(x,y,'r.',x,z)
title('Regular smoothing')
subplot(122), plot(x,y,'r.',x,zr)
title('Robust smoothing')

Effect of exaggeratedly high smoothness parameter: s is multiplied by 100.

z = smoothn(y,s*100);
plot(x,y,'r.',x,z,'k','LineWidth',2)

2-D smoothing

Create a noisy 2-D parametric surface and smooth it using smoothn.

[xp,yp] = deal(0:.02:1);
[x,y] = meshgrid(xp,yp);
f = exp(x+y) + sin((x-2*y)*3);
fn = f + randn(size(f))*0.5; % Adding white noise
fs = smoothn(fn); % DCT-based smoothing
subplot(121), surf(xp,yp,fn), zlim([0 8]), axis square
subplot(122), surf(xp,yp,fs), zlim([0 8]), axis square

2-D smoothing with missing data

n = 256;
y0 = peaks(n);
y = y0 + rand(size(y0))*2;
I = randperm(n^2);
y(I(1:n^2*0.5)) = NaN; % lose 1/2 of data
y(40:90,140:190) = NaN; % create a hole
z = smoothn(y); % smooth data
subplot(2,2,1:2), imagesc(y), axis equal off
title('Noisy corrupt data')
subplot(223), imagesc(z), axis equal off
title('Recovered data ...')
subplot(224), imagesc(y0), axis equal off
title('... compared with original data')

3-D smoothing

[x,y,z] = meshgrid(-2:.2:2,-2:.2:2,-2:.2:2);
xslice = [-0.8,1]; yslice = 2; zslice = [-2,0];
vn = x.*exp(-x.^2-y.^2-z.^2) + randn(size(x))*0.06;
subplot(121), slice(x,y,z,vn,xslice,yslice,zslice,'cubic')
title('Noisy data')
v = smoothn(vn);
subplot(122), slice(x,y,z,v,xslice,yslice,zslice,'cubic')
title('Smoothed data')

Cardioid

t = linspace(0,2*pi,1000);
x = 2*cos(t).*(1-cos(t)) + randn(size(t))*0.1;
y = 2*sin(t).*(1-cos(t)) + randn(size(t))*0.1;
z = smoothn(complex(x,y));
plot(x,y,'r.',real(z),imag(z),'k','linewidth',2)
axis equal tight

Reference

Garcia D, Robust smoothing of gridded data in one and higher dimensions with missing values. Computational Statistics & Data Analysis, 2009, doi:10.1016/j.csda.2009.09.020 (PDF).

About the author

Damien Garcia, Eng., Ph.D.
Assistant professor, Department of radiology
CRCHUM, University of Montreal Hospital
Montreal, QC, Canada
Damien.Garcia.REMOVE-THIS@BiomeCardio.com