Code covered by the BSD License  

Highlights from
Estimatenoise

5.0

5.0 | 11 ratings Rate this file 51 Downloads (last 30 days) File Size: 5.83 KB File ID: #16683

Estimatenoise

by

 

01 Oct 2007 (Updated )

Noise variance estimation from a signal vector or array

| Watch this File

File Information
Description

Some curve fitting or smoothing tools can benefit from knowledge of the noise variance to expect on your data. Kalman filters use this information, also some spline fitting tools. So I wrote a function to extract the noise variance from a signal vector. It also works on any specified dimension of an array.

A few examples of this code in use:

Simple linear data, with purely additive N(0,1) gaussian noise:

t = 0:10000;
x = t + randn(size(t));
mv = estimatenoise(x)
mv =
      1.0166

Gaussian noise added to a sine wave (Nominal variance = 0.01)
t = linspace(0,1,1000)';
x = sin(t*50) + randn(size(t))/10;
mv = estimatenoise(x)
mv =
    0.0096887

Pure gaussian noise, with a nominal variance of 9. (Note that var would have been a better estimator for this particular case...)

mv = estimatenoise(3*randn(2,3,1000),3)
mv =
       9.6584 8.2696 8.632
       9.2404 8.5346 9.7725

A piecewise constant function with multiple discontinuities. The true noise variance should be 0.01.

t = linspace(0,1,1000);
X = round(cos(t*6*pi)) + randn(size(t))/10;
plot(t,X)
var(X) % var will be wildly in error
ans =
     0.68256

estimatenoise(X)
ans =
      0.010882

Test if estimatenoise is able to recover the variance of a normally distributed random sample with unit variance. (Yes, it will be much slower than var.)

mean(estimatenoise(randn(1000,1000)))
ans =
       1.0002

Estimatenoise can now handle non-uniformly spaced series (by request.) In the next example,
the actual noise variance was 1.0 here. Perform the operation 1000 times, then look at the median variance estimate. How well did we do?

t = sort([randn(1,100) , randn(1,100)+5]);
X = repmat(sin(t*5)*100,1000,1) + ...
   randn(1000,length(t));

Estimatenoise is clearly wrong when the spacing is ignored.

median(estimatenoise(X,2))
ans =
      16.438

Supplying the sampling "times", we get quite a reasonable result.
median(estimatenoise(X,2,t))
ans =
      1.1307

Estimatenoise also works on data with replicates. In this example, each point will be replicated up to 3 times. The actual noise variance was again 1.0 here. I'll also compare the increase in times required for estimatenoise when t is supplied.

t = sort([0:1:100, 1:2:100, 1:4:100]);
X = repmat(sin(t/10)*100,1000,1)+ ...
  randn(1000,length(t));

Again, estimatenoise is clearly wrong when the non-uniform spacing is ignored.

tic,median(estimatenoise(X,2)),toc
ans =
      4.2056
Elapsed time is 2.690135 seconds.

Supplying the sampling "times", we again get quite a reasonable result. The time penalty is not quite 2x.

tic,median(estimatenoise(X,2,t)),toc
ans =
      1.0116
Elapsed time is 4.864486 seconds.

Acknowledgements

This file inspired Reverberation Time Calculator, Nth Octave Frequency Bands, Sound Power Directivity Analysis, Nth Octave Test Signal, Nth Oct Hand Arm & Ac Filter Tool Box, Continuous Sound And Vibration Analysis, and Noise Variance Estimation.

MATLAB release MATLAB 7.4 (R2007a)
Other requirements While written in version 7.4 (R2007a) this relatively simple code should run in older releases with no problems.
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (15)
23 Sep 2014 dang thanh

Please give me the paper/reference to explain this method! thanks

05 Jul 2014 Antoni J. CanĂ³s

Very nice. Thank you John!!

09 Sep 2011 lara presa

Thanks!

08 Sep 2011 Mark Shore

@ lara, John D'Errico recently announced he would be significantly reducing his input to MATLAB forums. This is a 1D tool, so you might want to consider using Damien Garcia's N-D noise estimation tool at http://www.mathworks.com/matlabcentral/fileexchange/25645

08 Sep 2011 lara presa

Can we estimate the noise in an image with this function?? what happen with sampling times t?? is necessary for a matrix? Thanks!!

06 Aug 2009 hbu 

thank you .

15 Dec 2008 Jody Klymak

^ Sorry, the first cumsum above should not be multiplied by 3!

14 Dec 2008 Jody Klymak

This is very nice.

Potential users should be aware that it assumes a spectral gap in the data. If you generate a red signal with noise:

>> x = cumsum(randn(1000,1000))*3 + randn(1000,1000);
>> mean(estimatenoise(x))
ans =

1.3021

>> x = cumsum(randn(1000,1000))*3 + randn(1000,1000);
>> mean(estimatenoise(x))
ans =

3.3083

you get a high biased estimate because there is still significant signal at high frequencies. As with any noise estimate, you need to understand the signal first!

I use this example because a large class of geophysical signals are red.

13 Dec 2008 Kenneth Eaton  
10 Dec 2007 John D'Errico

For those users who asked for the non-uniform version, its up there now. Sorry it took a while, I had a couple of aborted attempts before I got it right. The new code will use the old version on equally spaced series, but the unequally spaced option is not terribly slower in my tests.

03 Dec 2007 A B

I would like to see the extension to unevenly distributed data

10 Nov 2007 John de Mello

An extremely useful piece of code. The planned extension to unevenly spaced data sets would be fantastic.

28 Oct 2007 pule shule  
13 Oct 2007 Ihsan ulhaq

Yes well tested, verified, simple, mature and well documented code.
We should appreciate such job.

02 Oct 2007 G. A.M.

Very well documented. Very clear code. Nice examples. The algorithm has been well-tested (Although the author feels he could do more testing and verification and further improve this code, don't we all feel that way about code we write?) I think this code is more than mature enough for release here. I appreciate the author making it available.

Updates
18 Oct 2007

Update to handle complex inputs

Contact us