How to determine variance/std deviation of Gaussian noise from measured data

Question

If I have a signal measured in discrete time intervals of a process containing some signal and some Gaussian noise how do I go about measuring the variance of this noise [i.e. the $\sigma^2$ term in $\mathcal{N}(0, \sigma^2)$].

For example if I generate some pure Gaussian noise and plot the time trace, the histogram of counts and the power spectral density like so:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import scipy.signal

sigma = 0.1
sample_freq = 1e6 # 1 MS/s
dt = 1/sample_freq
r = np.random.normal(0, sigma, 1000000)
t = np.arange(0, 1000000*dt, dt)

# Create 2x2 sub plots
gs = gridspec.GridSpec(2, 2)

fig = plt.figure(figsize=[10, 7])
ax1 = fig.add_subplot(gs[0, 0]) # row 0, col 0
ax1.plot(t, r)

ax2 = fig.add_subplot(gs[0, 1]) # row 0, col 1
ax2.hist(r, bins=100, orientation="horizontal")
for tick in ax2.get_xticklabels():
    tick.set_rotation(270)

ax3 = fig.add_subplot(gs[1, :]) # row 1, span all columns
freqs, PSD = scipy.signal.welch(r, sample_freq, nperseg=1000)
PSD = PSD[freqs.argsort()]                                                  
freqs.sort()
ax3.semilogy(freqs, PSD)
ax3.set_ylabel("$S_{xx}$ ($V^2/Hz$)")
ax3.set_xlabel("$f$ (Hz)")

I could fit the equation for the probability density of a Gaussian [i.e. $f(x|\mu, \sigma^2) = \dfrac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{(\mu-x)^2}{2\sigma^2}}$] to the histogram to get the std deviation $\sigma$ of the noise.

However, in the case where we have some signal with some Gaussian noise we get something like that which I have below for a pure 100kHz sinusoidal signal:

rs = r + 1*np.sin(2*np.pi*100e3*t)

# Create 2x2 sub plots
gs = gridspec.GridSpec(2, 2)

fig = plt.figure(figsize=[10, 7])
ax1 = fig.add_subplot(gs[0, 0]) # row 0, col 0
ax1.plot(t, rs)

ax2 = fig.add_subplot(gs[0, 1]) # row 0, col 1
ax2.hist(rs, bins=100, orientation="horizontal")
for tick in ax2.get_xticklabels():
    tick.set_rotation(270)

ax3 = fig.add_subplot(gs[1, :]) # row 1, span all columns
freqs, PSD = scipy.signal.welch(rs, sample_freq, nperseg=1000)
PSD = PSD[freqs.argsort()]                                                  
freqs.sort()
ax3.semilogy(freqs, PSD)
ax3.set_ylabel("$S_{xx}$ ($V^2/Hz$)")
ax3.set_xlabel("$f$ (Hz)")

We can now no longer fit to the histogram without deconvolving the sinusoidal signal some how (which may not be so trivial for a more complex signal). However the noise level in the power spectral density remains the same, so there should be some way of getting the std deviation $\sigma$ of the Gaussian noise from this.

I understand that the equation for the power spectral density in discrete time is:

$$S_{xx}(\omega) = \dfrac{(\Delta t)^2}{T}|\sum_{n=1}^{N}x_n e^{-i\omega n}|^2$$

and the equation for the probability density of $x$ (where $x$ is drawn from $\mathcal{N}(0, \sigma^2)$) is

$$f(x|\mu, \sigma^2) = \dfrac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(\mu-x)^2}{2\sigma^2}}$$

But I am unsure how to derive the relationship between the level of the noise in the PSD and the std deviation $\sigma$ of the noise in order to extract the std deviation $\sigma$ in the noise on a measured signal.

Answer 1

What is typically done is to apply an orthogonal transform (or a close to tight frame) that

• preserves the Gaussianity,
• better separates the signal and the noise.

Fourier for simple sums of sines, wavelets for piecewise polynomial signals. Then, identify or parse some coefficient bins (bands, wavelet subbands), expecting that signal coefficients are sparse enough. Then, compute a robust scale estimator (like the Median Absolute Deviation, MAD) on the absolute values of the coefficients. And apply some factor to the result. This could work also with band-pass filters, derivative operators, etc.

You can find implementations in Matlab:

• wnoisest: Estimate noise of 1-D wavelet coefficients
• adapted for Ricean noise: Robust rician noise estimation for 3D MR

References:

• Nonlinear wavelet thresholding: A recursive method to determine the optimal denoising threshold