## Stochastic resonance¶

I'd like to adapt a figure from Zeng et al (2000) to illustrate a blog post about nonlinear audition. (Note, the post will be up at about 1300 UTC on Monday 9 June 2014.)

Zeng FG, Fu Q, Morse R (2000). Human hearing enhanced by noise. Brain Research 869, 251–255.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [2]:
x = np.arange(500)
y = 0.3 + np.sin(x/20.)/5
noisy = y + np.random.rand(500)/3 - (3./20)

In [3]:
plt.plot(x, y)
plt.plot(x, noisy, alpha=0.5)
plt.show()


## Build the figure¶

In [4]:
plt.figure(figsize=(12,4))

plt.subplot(121)
plt.plot(x, y, lw=2, color='k')
plt.text(10, 0.05, "signal we'd like to detect", color='k')
plt.hlines(0.6, 0, 500, color='b', linestyles='--')
plt.text(10, 0.62, 'detection threshold', color='b')
plt.ylabel('amplitude')
plt.xlabel('time (ms)')
plt.ylim(0, 0.8)
plt.xlim(0, 500)

plt.subplot(122)
plt.plot(x, noisy, color='g', alpha=0.8)
plt.plot(x, y, color='k', lw=2)
plt.hlines(0.6, 0, 500, color='b', linestyles='--')
plt.xlabel('time (ms)')
plt.ylim(0, 0.8)
plt.xlim(0, 500)

plt.show()


## 2D image example¶

Similarly, we can form a 2D signal below some threshold, infest it with noise, and try to recover it.

In [5]:
signal_level = 0.7
noise_level = 1.0
threshold_level = 1.1

In [6]:
im = np.zeros((64,64))
im[16:48,16:48] = signal_level

plt.imshow(im, interpolation="None")
plt.set_cmap('Greys')
plt.colorbar()
plt.show()


We'll define a function to threshold an array. Pass in an array (image) and a level to threshold at.

In [7]:
def threshold(im, value=0.5, mode='same'):
t = np.copy(im)
t[t<value] = 0
if mode == 'binary':
t[t>=value] = 1
return t


Now threshold so we can't see the square.

In [8]:
plt.imshow(threshold(im, value=threshold_level), interpolation="None")
plt.show()

In [9]:
noise = np.random.rand(64,64) * noise_level - noise_level/2
noised = im + noise
plt.imshow(noised, interpolation="None")
plt.colorbar()
plt.show()

In [10]:
noise = np.random.rand(64,64) * noise_level - noise_level/2
noised = im + noise
thresholded = threshold(noised, threshold_level)
plt.imshow(thresholded, interpolation="None", vmin=1)
plt.colorbar()
plt.show()

In [11]:
plt.figure(figsize=(18,4))

plt.subplot(141)
plt.title('Signal at {}'.format(signal_level))
plt.imshow(im, interpolation="None", cmap="Greys")
plt.ylim(0, 64)
plt.xlim(0, 64)
plt.hlines(32, 0, 64, color='r')

plt.subplot(142)
plt.title('Noise at {0}'.format(noise_level))
plt.imshow(noise, interpolation="None", cmap="Greys")
plt.ylim(0, 64)
plt.xlim(0, 64)
plt.hlines(32, 0, 64, color='r')

plt.subplot(143)
plt.title('Section across signal and noise', color='r')
plt.plot(im[32,:], lw=2, color='k')
plt.plot(noised[32,:], color='g')
plt.ylim(-0.1, 1.3)
plt.xlim(0, 64)
plt.hlines(1.1, 0, 64, color='b')
plt.text(3, 1.12, 'threshold', color='b')
plt.text(6, 0.63, 'signal', color='k')
plt.text(15, 0.90, 'noisy\nsignal', color='g', horizontalalignment='right')

plt.subplot(144)
plt.title('Detected signal')
plt.imshow(thresholded, interpolation="None", cmap="Greys", vmin=1)
plt.hlines(32, 0, 64, color='r')
plt.show()

plt.show()


We should probably add noise to the detected image. This is the thermal noise in the detection system, unlike the other noise, which was the noise we added to the 'input', so to speak.

In [12]:
rec_noise = np.random.rand(64,64) + 0.5
recorded = thresholded + rec_noise

plt.imshow(recorded, interpolation="None", cmap="Greys", vmin=1)
plt.show()


## Digging deeper¶

### Everything after this point is a work in progress and probably wrong.¶

What happens if we just add a ton of noise and threshold? Is that more realistic?

Also, let's switch to a different noise distribution. It seems that when I use a uniform distribution like numpy.random.rand I don't get the reduction in S:N that you expect when you add more and more noise. Instead, we'll use a logarithmic distribution. It's going to be a nasty hack to get reasonable standard deviation for it...

In [15]:
import scipy.stats

signal_level = 0.3
noise_level = 0.4
threshold_level = 1

im = np.zeros((64,64))
im[16:48,16:48] = signal_level

# Use a different distribution, with a nasty hack for the standard deviation
noise = np.random.lognormal(noise_level, (noise_level+1)/6, 64*64).reshape((64,64))
#noise = np.random.rand(64,64) * noise_level - noise_level/2
noised = im + noise

thresholded = threshold(noised, threshold_level)

plt.figure(figsize=(18,5))

plt.subplot(131)
plt.imshow(threshold(im, threshold_level), interpolation="None")
plt.title('S:N {0}'.format(scipy.stats.signaltonoise(threshold(im, threshold_level), axis=None)))
plt.colorbar(shrink=0.8)

plt.subplot(132)
plt.imshow(noised, interpolation="None")
plt.title('S:N {0}'.format(scipy.stats.signaltonoise(noised, axis=None)))
plt.colorbar(shrink=0.8)

plt.subplot(133)
plt.imshow(thresholded, interpolation="None")
plt.title('S:N {0}'.format(scipy.stats.signaltonoise(thresholded, axis=None)))
plt.colorbar(shrink=0.8)
plt.show()


We can track the S:N as we increase noise:

In [14]:
signal_level = 0.3
threshold_level = 1

im = np.zeros((64,64))
im[16:48,16:48] = signal_level

ston = []
r = np.arange(0,3.,0.05)

for n in r:
noise = np.random.lognormal(n, (n+1)/6, 64*64).reshape((64,64))
#noise = np.random.rand(64,64) * n - n/2
noised = im + noise
thresholded = threshold(noised, threshold_level)
ston.append(scipy.stats.signaltonoise(thresholded, axis=None))

plt.plot(r, ston)
plt.title('Max at {0}'.format(r[max(enumerate(ston),key=lambda x: x[1])[0]]))
plt.ylabel('signal:noise')
plt.xlabel('mean of noise distribution')
plt.show()


I think this looks about right. And we've established that the Goldilocks value is 0.55 for this type of noise.

But when I go back and explore the different levels of noise, it looks better for smaller values. Also, a perfect black square doesn't score well on this S:N measure. So I think I need to revisit this. Stay tuned...

In [14]: