In [3]:
from numpy import *
%pylab inline
Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['gamma', 'beta']
`%matplotlib` prevents importing * from pylab and numpy


In [33]:
from scipy.misc import factorial
from scipy.stats import poisson
#def p(x,N):
    #this naive definition will only work for small values on 32-bit install
#    return exp(-N)*(N**x)/factorial(x)
def p(x,N):
    return poisson.pmf(x, N)
t = arange(1,30,dtype=numpy.int64)
for N in [1,4,10,15]:
legend(loc='upper left')
<matplotlib.legend.Legend at 0x113d554d0>


Photons are generated by a light source independently and randomly with an average rate of $N$ per second. How many must be counted by a photodetector per second to determine the rate to within 1%? 1 ppm? What are the corresponding powers in watts for visible light?

Since the photons are random and independent, this is a Poisson process. We saw above that the relative standard deviation is given by $\frac{1}{\sqrt{N}}$, where $N$ is number of events. In our case, this represents detections per second. Thus, to determine the photons per second to within 1%, we must measure $.01^{-2} = 1\times10^{4}$ photons. To within 1 ppm, we must measure $(1\times10^{-6})^{-2} = 1\times10^{12}$.

A photon has energy given by $E = hc/\lambda$, where $\lambda$ is the wavelength. For visible light, say $\lambda = 500 nm$. Then $E = 4.0\times10^{-19} J$ or around $2 eV$.

In [10]:
h = 6.626e-34 #planck const
l = 500e-9 #m, wavelength
c = 3.0e8 #m/s, speed of light
ev = 1.6e-19 #J, electron volt
e_p = h*c/l
print e_p, e_p/ev
def power_print(N):
    print "The wattage of %.0e visible photons per second is %.2e watts."%(N,N*e_p)
3.9756e-19 2.48475
The wattage of 1e+04 visible photons per second is 3.98e-15 watts.
The wattage of 1e+12 visible photons per second is 3.98e-07 watts.


In [8]:
DF = 20e3 #20 kHz bandwidth
R = 10e3 #10 kOhm source impedance
T = 300 #Kelvin, room temp
snr = 20 #20 db desired snr wrt johnson noise
k = 1.38e-23
def v(df,t,r):
    return 10*sqrt(4*k*t*r*df)
def equiv_cap(v,t):
    #return capacitance in Farads for a capacitor exhibiting fluctuations of voltage v
    return k*t/v/v
t = linspace(T-50,T+50,100)
r = linspace(R-5e3,R+5e3,100)
df = linspace(DF-9e3,DF+9e3,100)
print "Minimum input voltage: %.2f uV"%(1e6*v(DF,T,R))
print "Johnson noise voltage: %.2f uV"%(1e6*v(DF,T,R)/10)
print "Equivalent capacitance: %.2f nF"%(1e9*equiv_cap(v(DF,T,R)/10,T))
subplot(1, 3, 1)
xlabel('temperature (K)')
ylabel('required input voltage ($\mu V$)')
subplot(1, 3, 2)
xlabel('resistance ($k\Omega$)')
subplot(1, 3, 3)
xlabel('amplifier bandwidth ($kHz$)')
Minimum input voltage: 18.20 uV
Johnson noise voltage: 1.82 uV
Equivalent capacitance: 1.25 nF
In [90]:
df=20e3 #bandwidth
q = 1.602e-19 #electron charge, Couloumb
def min_current(ratio, df):
    #return minimum current for ratio of current to rms of shot noise to be a given ratio
    #ratio: ratio of shot noise to current
    #df: bandwidth, Hz
    return 2*q*df/ratio/ratio
ratio = .01
print "Current for %.1f%% of rms shot noise to current is %.2e amps"%(100*ratio,min_current(ratio,df))
ri = logspace(-6,0,100)
xlabel("Ratio rms shot noise to current")
ylabel("Minimum current necessary (amps)")
Current for 1.0% of rms shot noise to current is 6.41e-11 amps


In [6]:
alpha = .05
beta = .05
gamma = alpha+beta
c1 = 1./gamma
c2 = alpha/gamma
def p(t):
    return asarray(c1*array([beta,alpha]) + c2*exp(-gamma*t[...,None])*array([1,-1]))
def sigma(t):
    mean = alpha/gamma - c2*exp(-gamma*t)
    return sqrt(mean-mean**2)
def R(tau):
    autocovariance = (alpha**2-c2*alpha)/gamma**2 + (c2**2/(2*gamma)+c2*alpha/gamma/gamma)*exp(-gamma*tau)
    return autocovariance#/sigma(tau)
t = linspace(0,20,200)
ax1 = plt.gca()#plt.subplot(211)
ax1.legend(loc='lower right')

#ax2 = plt.subplot(212, sharex=ax1)
#ax1.set_title(r'probabilaties, $\alpha=%.2f,\hspace{.5}\beta=%.2f$'%(alpha,beta))
(0, 1)

For the autocorrelation, we can use these functions to compute the expectation of $\lt x(t) x(t-\tau) \gt$. The expression is only nonzero when both $x(t)$ and $x(t-\tau)$ are nonzero. So, given that $x(t-\tau)$ equals 1, what is probability that $x(t)$ also equals 1? This is exactly the expression we derived above for $p_1$ : $\frac{1}{\alpha+\beta}(\alpha + \beta e^{-(\alpha+\beta)\tau} ) $. This has the right behavior at $\tau=0$ and at steady state.

Since we're taking a limit for the expectation, we can use the steady state probability that the system is in state 1: $\frac{\alpha}{\alpha+\beta}$ as the prior.

Thus, the autocorrelation function is $\lt x(t) x(t-\tau) \gt = \frac{\alpha}{(\alpha+\beta)^2}(\alpha + \beta e^{-(\alpha+\beta)\tau} ) $.

To compute the power spectrum, we need to take a Fourier transform of this autocorrelation (Wiener-Khinchin theroem). Let's just consider a Fourier transform of the form $A+Be^{-\gamma t}$. Fourier is linear, so we can just consider the transorm of a constant and an exponential. The transform of a constant is a delta function. The transform of an exponential can be computed as follows:

$F[e^{-\gamma t}] = \int_0^{\infty} (\cos(2\pi \gamma t) + i \sin(2\pi\gamma t)) e^{-2\pi\gamma t} dt + \int_0^{\infty} (\cos(2\pi \gamma t) - i \sin(2\pi\gamma t)) e^{-2\pi\gamma t} dt$

$F[e^{-\gamma t}] = 2\int_0^{\infty} \cos(2\pi\gamma t)e^{-2\pi\gamma t} dt$

Punting, wolfram says this can be computed by integration by parts:

$F[e^{-\gamma t}](f) = \frac{1}{\pi} \frac{\gamma}{\gamma^2 + f^2}$

Which is a Lorentzian where $\tau = \gamma/2 = \frac{\alpha+\beta}{2}$ in our case.

In [ ]:
In [ ]: