"""
shot_noise.py
OVERVIEW
This Python program simulates shot noise such as that associated with the flow
of electrons across a reverse-biased PN junction. Shot noise is a type of
electronic noise originating from the discrete nature of electric charge.
Our model is as follows:
- Electrons arrive at and begin crossing the junction at random times; these
times can be modeled as a Poisson process with specified rate lambda. The
Poisson assumption implies that the number of electrons arriving during any time
interval is independent of the number arriving during any non-overlapping time
interval.
- The time T_transit for an electron to transit the junction is modeled as a
constant. In actuality, the transit time varies, but for purposes of an initial
investigation, this is a reasonable assumption. [Allowing the transit time to
have a distribution with a tail causes the autocorrelation to decay smoothly to
zero, and eliminates the nulls in the power spectral density.]
The current in units of Amperes is given by the following equation:
I(t)= n(t) * q_e / T_transit
In the above, n(t) is the number of electrons in transit (moving through the
junction) at time t, q_e is the charge of one electron, and T_transit is the
time for an electron to pass through the junction. Thus, the current is
proportional to the number of electrons in transit.
OBJECTIVES
The objectives of the simulation are to verify the following:
(1) The distribution of the number of electrons in transit should be
approximately normal (Gaussian).
(2) The autocorrelation should have a triangular shape, diminishing to zero when
the lag reaches T_transit.
(3) The power spectral density should be nearly white up to frequencies on the
order of 0.25 / T_transit.
IMPLEMENTATION APPROACH
Because electrons do not interact with one another in our model, we can generate
all of the arrival and departure times at the outset and then calculate the
number of electrons in transit at each event time. The use of a DES framework
such as SimPy would be overkill for this model. In detail, the following
sequence of calculations is performed:
(1) Generate the total number of electrons that arrive over the course of the
simulation and store this value as `N_arrivals`. This number is a Poisson
random variable with mean lambda * T_sim.
(2) Generate the arrival times. Given the number of arrivals and the total
simulation time, the arrival times are i.i.d. uniform on [0, T_sim].
(3) Add the transit time `T_transit` to each arrival time to obtain the
corresponding departure time. (This is done as a single array operation to
avoid slow, Python-level looping).
(4) Combine the lists of arrival and departure times, add the time zero to the
combined list, and sort the result to obtain a list of all event times of
interest. These are stored in the array `times`.
(5) Remove elements of `times` that are greater than or equal to `T_sim`.
(6) For each time in `times`, calculate the number of electrons in transit,
i.e., inside the junction. Note that this is the number of arrivals minus the
number of departures.
CONCLUSIONS
The histogram confirms that the distribution of electrons in transit is
approximately normal. (It becomes closer to normal as the mean number of
electrons in transit increases).
The autocorrelation graph has the expected triangular shape, except for some
noise that could be reduced by increasing T_sim.
The PSD is not as flat as expected at low frequencies. Probable causes are
either insufficient data or limitations of the PSD estimation algorithm.
REFERENCES
1. For general background on thermal noise and shot noise:
http://www.fis.unipr.it/~gigi/dida/strumentazione/harvard_noise.pdf
2. For background on the Poisson process:
http://en.wikipedia.org/wiki/Poisson_process
3. Bandwidth of a transistor depends on junction transit time, which in turn
depends on the materials, geometry, and temperature (higher temperature reduces
the junction transit time). See, for example, the following:
http://ecee.colorado.edu/~bart/book/book/chapter5/ch5_5.htm
AUTHOR
Dr. Phillip M. Feldman
LAST UPDATE
1 Jan, 2015
"""
# Section 1: Import functions and classes from modules.
import matplotlib, numpy, pdb
from numpy import array, cos, pi, sum
from numpy.random import rand, RandomState
from console_input import get_float, get_int, make_check
from matplotlib import pyplot
from scipy.stats import norm
# Section 2: Define utility functions. Most of these are window functions for
# use with `matplotlib.mlab.psd`.
def blackman(x):
return numpy.blackman(len(x)) * x
def hamming(x):
return numpy.hamming(len(x)) * x
def hanning(x):
return numpy.hanning(len(x)) * x
def nuttall(x):
N= len(x)
n= numpy.arange(N)
return 0.355768 - 0.487396*cos(2*pi*n/(N-1)) - 0.144232*cos(4*pi*n/(N-1)) \
- 0.012604*cos(6*pi*n/(N-1))
def autocorr(x):
"""
This function calculates the autocorrelation of a sequence.
"""
result= numpy.correlate(x, x, mode='full')
result= result[result.size/2:]
return result / result[0]
# Section 3: Get user inputs and initialize random number generator.
# In the following statement, we can't use `lambda` as a variable name because
# this is a Python reserved word. Appropriate choice of `get_float` or
# `get_int`, together with a suitable `check` function, prevents the user from
# entering an illegal value.
lam = get_float("\nInput lambda= arrival rate of electrons "
"(default=30.0): ", check=make_check('x>0'), default= 30.0)
T_transit= get_float("\nInput T_transit= time for an electron to transit the "
"junction (default=1.0): ", check=make_check('x>0'), default= 1.0)
T_sim = get_float("\nInput T_sim= total simulation time (default=3000.0): ",
check=make_check('x>0'), default= 3000.0)
# Create random number generator stream and initialize it. (In this case we
# need only a single stream, but some simulation applications require multiple
# streams).
RNG= RandomState()
# Section 4: Generate arrival and departure times. (The calculations performed
# here correspond to steps (1)-(5) from the list in the header document).
# Draw random number of arrivals from Poisson distribution with parameter
# lambda * T_sim:
N_arrivals= RNG.poisson(lam*T_sim)
# Generate the arrival times. Given the number of arrivals and the total
# simulation time, the arrival times are i.i.d. uniform on [0, T_sim]:
arrival_times= numpy.sort(rand(N_arrivals)*T_sim)
departure_times= arrival_times + T_transit
# Create array `times` that includes all arrival times, all departure times, and
# zero, with any duplicates removed:
times= numpy.unique( [0.0] + list(arrival_times) + list(departure_times) )
# Remove times greater than or equal to `T_sim`:
times= times[times < T_sim]
# Section 5: For each time in `times`, calculate the number of electrons in
# transit, i.e., inside the junction. Note that this is the number of arrivals
# minus the number of departures. (This calculation corresponds to step (6)
# from the list in the header document).
# The following calculation method works but is extremely wasteful of memory:
# in_transit= sum( arrival_times [:,None] <= times[None,:], axis=0 ) \
# - sum( departure_times[:,None] <= times[None,:], axis=0 )
# The following memory-efficient calculations uses `numpy.searchsorted`:
in_transit= numpy.searchsorted(arrival_times , times, 'right') \
- numpy.searchsorted(departure_times, times, 'right')
# Section 6: Plot the number of electrons in transit versus time.
print("Generating figure #1 ...")
# Open a new figure window, setting the face color to white:
fig= pyplot.figure(figsize=(10,6), dpi=120, facecolor=[1,1,1])
# Create axes object for single sub-plot that fills the entire window:
axes= fig.add_subplot(111)
# Generate step plot.
pyplot.step(times, in_transit, linewidth=1, where='post')
# Label the axes:
axes.set_xlabel('Time', size=16)
axes.set_ylabel('Number of Electrons in Transit', size=16)
# Turn on the grid:
pyplot.grid(True)
# Save figure to disk file:
pyplot.savefig('in_transit_vs_time.png', bbox_inches=0)
# Section 7: Delete the first 1 percent of the data to remove the initial
# transient.
i= int( 0.01 * len(in_transit) )
in_transit= in_transit[i:]
# Section 8: Plot a histogram of the number of electrons in the junction, with
# an overlaid normal curve (suitably scaled).
print("Generating figure #2 ...")
# Open a new figure window, setting the face color to white:
fig= pyplot.figure(figsize=(9,6), dpi=120, facecolor=[1,1,1])
# Create axes object for single sub-plot that fills the entire window:
axes= fig.add_subplot(111)
x_min, x_max= 10, 55
counts, bins= numpy.histogram(in_transit, bins=range(x_min, x_max+1))
pyplot.hist(in_transit, bins=bins)
pyplot.hold(True)
mu = numpy.mean(in_transit, dtype=numpy.float64)
sigma= numpy.std (in_transit, dtype=numpy.float64)
x= bins[:-1] + 0.5
y= counts.sum() * norm.pdf(x, loc=mu, scale=sigma)
pyplot.plot(x, y, linewidth=2, color='red')
pyplot.xlim([x_min, x_max])
axes.set_xlabel('Number of Electrons in Transit', size=16)
axes.set_ylabel('Counts', size=16)
# Save figure to disk file:
pyplot.savefig('in_transit_histogram.png', bbox_inches=0)
# Section 9: Using the original arrival and departure time data, generate a new
# array storing the number of electrons in transit for uniformly-spaced times.
# (To generate a power spectral density using conventional methods, we must have
# uniformly-spaced times).
# Oversample to minimize aliasing artifacts:
T_sample= 1.0 / (4.0 * lam)
times= numpy.linspace(0.0, T_sim, T_sim/T_sample)
# For each time in `t`, calculate the number of electrons in transit, i.e.,
# inside the junction. Note that this is the number of arrivals minus the
# number of departures.
in_transit= numpy.searchsorted(arrival_times , times, 'right') \
- numpy.searchsorted(departure_times, times, 'right')
# Delete the first 1 percent of the data to remove the initial transient:
i= int( 0.01 * len(in_transit) )
in_transit= in_transit[i:]
# Truncate array of times to the same length as `in_transit`:
times = times[:len(in_transit)]
# Remove the mean:
in_transit= in_transit - in_transit.mean()
print("Length of sequence to be used for autocorrelation and PSD: %d"
% len(in_transit))
# Section 10: Plot the sample autocorrelation of the number of electrons in the
# junction.
print("Generating figure #3 (this may take up to 60 seconds) ...")
# Open a new figure window, setting the face color to white:
fig= pyplot.figure(figsize=(10,6), dpi=120, facecolor=[1,1,1])
# Create axes object for single sub-plot that fills the entire window:
axes= fig.add_subplot(111)
pyplot.plot(times, autocorr(in_transit), linewidth=2)
pyplot.xlim([ 0.0, 10.0*T_transit])
pyplot.ylim([-0.1, 1.0])
axes.set_xlabel('Time Difference', size=16)
axes.set_ylabel('Autocorrelation', size=16)
axes.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%.3g'))
# Set x-axis tick interval to unity:
axes.xaxis.set_major_locator (matplotlib.ticker.MultipleLocator(1.0))
pyplot.grid(True)
# Save figure to disk file:
pyplot.savefig('autocorrelation.png', bbox_inches=0)
# Section 11: Plot an estimate of the power spectral density of the number of
# electrons in the junction.
print("Generating figure #4 ...")
# Open a new figure window, setting the face color to white:
fig= pyplot.figure(figsize=(10,6), dpi=120, facecolor=[1,1,1])
# Create axes object for single sub-plot that fills the entire window:
axes= fig.add_subplot(111)
windows= {'Nuttall':nuttall, 'Hanning':hanning}
for NFFT in [2048, 8192]:
for key in windows:
# Calculate PSD and normalize peak value to unity:
P_xx, freqs= matplotlib.mlab.psd(in_transit, NFFT=NFFT,
Fs=1.0/T_sample, noverlap=3*NFFT/4, window=windows[key])
P_xx= P_xx / P_xx.max()
f= freqs[freqs <= 2.0]
pyplot.plot(f, P_xx[:len(f)], linewidth=2, label='%d/%s' % (NFFT, key))
pyplot.hold(True)
pyplot.ylim([1.e-4, 1.01])
axes.set_xlabel('Frequency', size=16)
axes.set_ylabel('Power Spectral Density', size=16)
axes.set_yscale('log')
axes.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%.3g'))
pyplot.grid(True)
pyplot.legend(loc='lower left')
# Save figure to disk file:
pyplot.savefig('power_spectral_density.png', bbox_inches=0)
pyplot.show()