Noises and Signals

ECE 6323

Only for students of ECE 6323/5358/6358                H. Q. Le Copyright

1. Noise basic concept and Power Spectral Density function

1.1 Introduction

1.1.1 Example 1

We measure some quantity as a function of time or space, or whatever. Our measurement is a set of numbers like this: (it has 50 elements):
{7,6,12,9,15,13,9,11,9,10,10,6,9,11,11,8,10,11,10,3,17,12,12,10,13,11,15,7,14,
  13,10,10,9,9,18,10,6,12,9,7,10,6,13,7,10,9,6,11,12,12};
We can plot our measurement like this:

Manipulate[
If[TrueQ[GENERATE],
x1=RandomVariate[NormalDistribution[10.,3],100] ;
x2=RandomVariate[NormalDistribution[10.,0.6],100] ;
x1plot=ListPlot[x1, Joined-> True, PlotStyle-> {Blue,Thickness[0.01]}
,PlotRange->{0.,20}, ImageSize-> 350] ;
x2plot=ListPlot[x2, Joined->True, PlotStyle->{Red,Thickness[0.01]}
,PlotRange->{0.,20}, ImageSize-> 350];
GENERATE=False,] ;

If[TrueQ[STAT==1],
{mean,var,SD} ={Mean[#],Variance[#],StandardDeviation[#]}& @ x;
xstatg=Grid[
MapThread[Prepend,{
Join[{{"Mean","Variance","StandardDeviation"}}
,{{mean,var,SD}
}]
,{" ","x"}}]
(* style  *)
,BaseStyle-> {20,Black,FontFamily-> "Calibri"}
,Background->{{Hue[0.5,0.4,1]},{Hue[0.12,0.71,1]}}
,Frame-> All,Spacings->{1,0.5}];
statgraphLine=Graphics[{Hue[0.32,1.,1.,0.15],Rectangle[{1,mean-SD},{Length[x],mean+SD}]
,Black,Opacity[1],Thickness[0.008],Line[{{1,mean},{Length[x],mean}}]
,Hue[0.85,1.,1.] ,Dashed,Line[{{1,mean-SD},{Length[x],mean-SD}}]
,Line[{{1,mean+SD},{Length[x],mean+SD}}]}];
If[TrueQ[HIST],
{xmin,xmax}={Min[x],Max[x]};
histy=BinCounts[x,{xmin,xmax,Abs[xmax-xmin]10./Length[x]}];
histmax=1.3*Max[histy] ;
statgraphHist=Graphics[{Hue[0.32,1.,1.,0.15],Rectangle[{mean-SD,0.},{mean+SD,histmax}]
,Black,Opacity[1],Thickness[0.008],Line[{{mean,0.},{mean,histmax}}]
,Hue[0.85,1.,1.] ,Dashed,Line[{{mean-SD,0.},{mean-SD,histmax}}]
,Line[{{mean+SD,0.},{mean+SD,histmax}}]}], ];

statval=Grid[
MapThread[Prepend,{
Join[{{"Mean","Variance","StandardDeviation"}}
,{{mean,var,SD}

}]
,{" ","x"}}]
(* style  *)
,BaseStyle-> {20,Black,FontFamily-> "Calibri"}
,Background->{{Hue[0.5,0.4,1]},{Hue[0.12,0.71,1]}}
,Frame-> All,Spacings->{1,0.5}] ;


STAT=Mod[STAT+1,3],
If[TrueQ[STAT==0],statgraphLine= statgraphLine0; statgraphHist=statgraphHist0 ;,] ] ;

{Show[xplot,statgraphLine]
,If[TrueQ[HIST],
{xmin,xmax}={Min[x],Max[x]};
histy=BinCounts[x,{xmin,xmax,Abs[xmax-xmin]10./Length[x]}];
Show[Histogram[x,PlotRange->{All,All},  Frame-> True, GridLines-> Automatic, ChartStyle->{Hue[0.55,1.,1.,0.5]}, ImageSize-> 350]

(* Histogram[x, {0,20,1},PlotRange->{{0,20},All},  Frame-> True, GridLines-> Automatic, ChartStyle->{Hue[0.55,1.,1.,0.5]}, ImageSize-> 350] *)
,statgraphHist],]
,If[STAT>0 && TrueQ[GRIDVAL],statval,]}
,Column[{Grid[{{Button["x1",{If[NumberQ[X1],X1=Mod[X1+1,2],X1=1]},BaseStyle-> {24, Bold,GrayLevel[Dynamic[0.7(1-X1)]],FontFamily-> "Calibri"}
,Background-> Hue[0.03,0.5,0.5(1+Dynamic[X1])]
]
,Button["x2",{If[NumberQ[X2],X2=Mod[X2+1,2],X2=1]},BaseStyle-> {24, Bold,GrayLevel[Dynamic[0.7(1-X2)]],FontFamily-> "Calibri"}
,Background-> Hue[0.25,0.5,0.5(1+Dynamic[X2])]
]
,Button["y1",{If[NumberQ[Y1],Y1=Mod[Y1+1,2],Y1=1]},BaseStyle-> {24, Bold,GrayLevel[Dynamic[0.7(1-Y1)]],FontFamily-> "Calibri"}
,Background-> Hue[0.55,0.5,0.5(1+Dynamic[Y1])]
]
,Button["y2",{If[NumberQ[Y2],Y2=Mod[Y2+1,2],Y2=1]},BaseStyle-> {24, Bold,GrayLevel[Dynamic[0.7(1-Y2)]],FontFamily-> "Calibri"}
,Background-> Hue[0.75,0.5,0.5(1+Dynamic[Y2])]
]}}]
,Grid[{{Button[,{GENERATE=True}
(*,BaseStyle-> {18, Bold,FontFamily-> "Calibri"}
,Background-> Hue[0.05,1,1] *)]

,Button["Histogram",{If[TrueQ[HIST],HIST=False,HIST=True]},BaseStyle-> {16, Bold,FontFamily-> "Calibri"}
,Background-> Hue[0.8,0.35,1]]
,Button["Statistics",{If[NumberQ[STAT],STAT=Mod[STAT+1,3],STAT=1]},BaseStyle-> {16, Bold,FontFamily-> "Calibri"}
,Background-> Hue[0.17,1,1]]
,Button["Numerical",{If[TrueQ[GRIDVAL],GRIDVAL=False,GRIDVAL=True]},BaseStyle-> {16, Bold,FontFamily-> "Calibri"}
,Background-> Hue[0.27,1,1]]}}]
}]
,Initialization:>{
statgraphLine0=Graphics[Point[{0.,0}]] ;
statgraphLine=statgraphLine0 ;
statgraphHist0=Graphics[Point[{0.,0}]] ;
statgraphHist=statgraphHist0 ;
X1=1; X2=0; Y1=0;Y2=0;
}]

The signal is not constant, but fluctuates a bit. What is the meaning of this fluctuation? is it meaningful? or meaningless? Does it contain information or not? There is no way we can answer this question beyond a shadow of doubt. But if we don't find this variation with any useful information to us, we call it noise.

So, we define the concept of noise, but how do we describe it? Suppose we do the measurement again and this time with get a new result. Let's compare the previous signal with this new result:

xnew=RandomVariate[NormalDistribution[10.,0.6],100] ;
xnewplot=ListPlot[xnew, Joined->True, PlotStyle->{RGBColor[1,0,0],Thickness[0.01]}
,PlotRange->{0.,20}, DisplayFunction->Identity] ;
Show[xplot,xnewplot]

We would intuitively say that the red is less noisy than the blue. But how much is less noisy? how do we quantify the "noisiness"?

1.2 Statistical perspective

Suppose that we KNOW THAT THE QUANTITY WE MEASURE IS SUPPOSED TO BE CONSTANT, then the fluctuation we see is definitely by definition, noise. We can plot the histogram like this

Row[{Histogram[{x,xnew}, {0,20,1},PlotRange->{{0,20},All},  Frame-> True, GridLines-> Automatic, ChartStyle->{Blue,Red}, ImageSize-> 400]
,Show[xplot,xnewplot, ImageSize-> 400]}]

They are so different. We can estimate the mean and variance

Grid[
MapThread[Prepend,{
Join[{{"Mean","Variance","StandardDeviation"}}
,{{Mean[x],Variance[x],StandardDeviation[x]}
,{Mean[xnew],Variance[xnew],StandardDeviation[xnew]}
}]
,{" ","x","xnew"}}]
(* style  *)
,BaseStyle-> {20,Black,FontFamily-> "Calibri"}
,Background->{{Hue[0.5,0.4,1]},{Hue[0.12,0.71,1]}}
,Frame-> All,Spacings->{1,0.5}]

Clearly, both signals have very similar mean, but one has a much smaller variance than the other. So, in this case, variance is a measure of the noise.

Can we always use variance as a measure of the noise?

Example 2

Suppose we have a new signal that looks like this:
y1={10.2,12.7,15.4,17.6,19.1,20.4,19.3,18.6,17.,14.2,11.5,8.65,6.05,3.52,1.52,
  0.46,0.294,1.03,2.19,4.67,7.02,9.62,13.4,15.2,18.1,19.1,20.,20.,18.1,16.6,
  14.,10.7,7.97,5.35,2.61,1.16,-0.129,0.2,0.684,2.48,4.43,8.04,11.,13.6,15.8,
  18.,19.4,20.1,19.7,18.} ;
  and
  y2={7.,6.3,12.6,9.78,15.9,14.,9.97,11.9,9.68,10.4,10.1,5.84,8.56,10.3,10.1,7.02,
  9.,10.1,9.23,2.45,16.7,12.,12.3,10.6,13.8,11.9,16.,7.97,14.9,13.7,10.4,10.1,
  8.83,8.54,17.3,9.12,5.02,11.,8.08,6.24,9.46,5.74,13.,7.33,10.6,9.8,6.94,12.,
  13.,12.8};
  we can find their variances:

y1=10*Table[Sin[2 Pi (i-1)/50.],{i,100}]+RandomVariate[NormalDistribution[10.,0.5],100] ;
y2= RandomVariate[NormalDistribution[10.,3],100] ;

Grid[
MapThread[Prepend,{
Join[{{"Mean","Variance","StandardDeviation"}}
,Table[{Mean[#],Variance[#],StandardDeviation[#]}& @ ys
,{ys,{y1,y2}}]
]
,{" ","y1","y2"}}]
(* style  *)
,BaseStyle-> {20,Black,FontFamily-> "Calibri"}
,Background->{{Hue[0.5,0.4,1]},{Hue[0.12,0.71,1]}}
,Frame-> All,Spacings->{1,0.5}]

So should we say that y1 is noisier than y2?

Row[{Histogram[{y1,y2}, {0,20,1},PlotRange->{{0,20},All},  Frame-> True, GridLines-> Automatic, ChartStyle->{Blue,Red}, ImageSize-> 400]
,Show[xplot,xnewplot, ImageSize-> 400]}]

yplot=ListPlot[{y1,y2},Joined->True, PlotStyle->{Thickness[0.01],Red,Blue}]

In the above, y1 is blue (variance= 50) and y2 is red variance is 9.8. Their means are:

Mean[y1]
Mean[y2]

Again, similar mean. Is it correct to say that y1 is noisier than y2?
This depends! If both y1 and y2 are supposed to be constant, we say that y1 fluctutates more and thus, noisier than y2. On the other hand, it is clearly that the fluctuation of y1 is somewhat predictable, and if we remove this "predictable" fluctuation, it may be less noisy. Thus, this is a key concept: noise BY DEFINITION CAN NOT BE DETERMINISTIC: in other words, noise has to be defined as the unpredictable, indeterministic part of the signal. One way to see in this example is to Fourier transform the signal.

fy1=ListPlot[Take[10 Log[10,((Abs[Fourier[y1-Mean[y1]]])^2)],{1,Ceiling[Length[y1]/2]}],Joined->True,PlotStyle->{RGBColor[0,0,1],Thickness[0.01`]},PlotRange->{-20,30},Frame->True,ImageSize->{300,300},AspectRatio->1,GridLines->Automatic,DisplayFunction->Identity];
fy2=ListPlot[Take[10 Log[10.`,((Abs[Fourier[y2-Mean[y2]]])^2)],{1,Ceiling[Length[y1]/2]}],Joined->True,PlotStyle->{RGBColor[1,0,0],Thickness[0.01`]},DisplayFunction->Identity];
Show[fy1,fy2,PlotRange-> All]

What we plot here is the Fourier transform square, on Log scale (unit dB) of the signal. It is clear here, that the blue (y1) has  a periodic signal at bin ~ 2 that is much stronger than that of the red, and the y1 F-transform signal at higher frequency is actually less than that of y2. In other words, the blue (y1) has a large variance, but seems less noisy than the red (y2).
We see how useful Fourier transform is.
Let's do the follow: we define a new signal y1new that is from y1 but subtracted by a known function that we guess to be a sine function:

y1new=y1-Table[10*Sin[x*0.3],{x,0,49,1}]

Now we compare y1new and y2:

y1newplot=ListPlot[y1new,Joined->True,PlotStyle->{RGBColor[0,0,1],Thickness[0.01`]},PlotRange->{0.`,20},DisplayFunction->Identity];
y2plot=ListPlot[y2,Joined->True,PlotStyle->{RGBColor[1,0,0],Thickness[0.01`]},PlotRange->{0.`,20},DisplayFunction->Identity];
Show[y1newplot,y2plot,DisplayFunction->$DisplayFunction]

From this, we have to say that y1 is less noise than y2! because eventhough it has a large fluctuation, we deterministically know exactly one component of that large fluctuation, (the sine function) and thus, it can be removed.

Example 3

Now we have a fairly complicated signal that looks like this:

x1=RandomReal[NormalDistribution[0,1],2000];
sig1=x1+Table[0.5` Sin[x 0.7`]+1.3` Sin[x 0.13`],{x,0,1999,1}];
ListPlot[sig1,Joined->True, AspectRatio-> 0.2]

How do we analyze this? well, we can Fourier transform again:

ListPlot[Take [10*Log[10, Abs[Fourier[sig1-Mean[sig1] ]]^2],{2,1000} ]
		    ,Joined->True, 
	PlotStyle->{RGBColor[0,0,1],Thickness[0.002]}
		 , PlotRange->{-20,30}
		, Frame->True, ImageSize->{650,300}, GridLines->Automatic]

What do we see? is there a signal in there? how much is the noise. Intuitively, we see 2 big peaks that stand out above the rest. We would call those peak signals, and the rest, noise. But what are we doing here, plotting the square of the Fourier transform of the signal? what is the meaning of this? Why are we doing this?

1.3 Power spectral density function

1.3.1 Definition

We define the power-spectral-density (PSD) function of a signal s(t) as:  
           N(f)= (lim_(T→∞)) (1/T)|(∫_0^T)s(t) (e^(i 2π f t)) dt  ((|)^2)

This is an important definition. To see what it does (or the meaning of it), let the signal be:  
                               s(t) = A ((e )^(-i 2 π t (Subscript[f, 0]) ));
then 
 N(f)= (lim_(T→∞)) (1/T)|(∫_0^T)s(t) (e^(i 2π f t)) dt ((|)^2)= (lim_(T→∞)) (1/T)|(∫_0^T)A ((e )^(-i 2 π (Subscript[f, 0]))) (e^(i 2π f t)) dt((|)^2) 
                                                 (Subscript[=  lim, T→∞]) ((|A(|^2))/T)|(∫_0^T) (e^(i 2π (f -(Subscript[f, 0]) ) t)) dt((|)^2)

FullSimplify[(∫_0^T) (E^(I 2π (f -(Subscript[f, 0]) ) t)) \[DifferentialD]t]

Now, we take the magnitude squared:

((FullSimplify[Expand[(-1+(E^(2 I π T (f-(Subscript[f, 0])))))(-1+(E^(-2 I π T (f-(Subscript[f, 0])))))]])/(4 (π^2) (((f-(Subscript[f, 0])))^2)))

(∫_(-∞)^∞)(1/T)(((((sin(π T f))/(π  f))))^2)\[DifferentialD]f

Therefore:
                    N(f)=  (lim_(T→∞)) ((|A(|^2))/T)((|(∫_0^T) (e^(i 2π (f -(Subscript[f, 0]) ) t)) dt|)^2)
But:            (lim_(T→∞)) (1/T)((|((sin(π T (f-(Subscript[f, 0]))))/(π  (f-(Subscript[f, 0]))))|)^2)  =   δ(f-(Subscript[f, 0])) 
Hence:           N(f)=(Subscript[lim , T→∞]) ((|A(|^2))/T)((|((sin(π T (f-(Subscript[f, 0]))))/(π  (f-(Subscript[f, 0]))))|)^2)= |A(|^2) δ(f-(Subscript[f, 0]))  
In other words, the PSD function of this signal tells us that it is zero for all frequency except for (Subscript[f, 0]) where it is a delta function with an amplitude factor (((|A|))^2).

That's why we call it PSD: it is a measure of the density of the spectral content of the signal as a function of frequency. 
The total power is finite: 
                   (∫_(-∞)^∞)N(f)df = |A(|^2) ∫δ(f-(Subscript[f, 0])) df  =  |A(|^2)  
(the power density is infinite at f=(Subscript[f, 0]) because the signal has purely (Subscript[f, 0]) frequency)

1.3.2 Practical calculation of PSD

If we have a finite series of signal measurement (or sampled points): (Subscript[{(Subscript[s, k])}, k=1]^M), how do we find its PSD?
Let say the series (Subscript[{(Subscript[s, k])}, k=1]^M) are measured with regular interval (either time, space, or whatever variable we used), then the series (Subscript[{(Subscript[s, k])}, k=1]^M) actually represent pairs of number:            (Subscript[{((Subscript[t, k]) , (Subscript[s, k]))  }, k=1]^M)  ≡(Subscript[{( (k-1)Δt, (Subscript[s, k]))  }, k=1]^(M  ))
where Δt is the sampling interval. We can take the Fourier transform of this numerical series by approx. the integral as a summation:
(∫_0^(M Δt))s(t)(e^(i 2π f t))dt  ≈ (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π f (k-1) Δt))  Δt
                                        = Δt (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π f (k-1) Δt))  
But if we discretize the variable f also, as a series (p-1) Δf for integer p, and choosing  Δf =1/(M Δt), we have the following expression:
(∫_0^(M Δt))s(t)(e^(i 2π  (Subscript[f, p]) t))dt  ≈ Δt (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π (Subscript[f, p]) (k-1) Δt))
                                                      = Δt  (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π (p-1) (k-1) Δt/(M Δt)))  
                                                      = Δt  (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π (p-1) (k-1) /M))
But we recall that in fast Fourier transform numerical technique, the Fourier transform (Subscript[b, s]) of a list (Subscript[a, r]) of length n is defined to be (1/(Sqrt[n]))(∑_(r=1)^n)(Subscript[a, r])(e^(2πi(r-1)(s-1)/n))
Therefore, the above expression can be written as:
(∫_0^(M Δt))s(t)(e^(i 2π  (Subscript[f, p]) t))dt  ≈ Δt  (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π (p-1) (k-1) /M))
                                          = Δt  (Sqrt[M]) FFT((Subscript[{(Subscript[s, k])}, k=1]^M);  p)

This is the basis for numerical calculation. We thus obtain the PSD as: 
 N((Subscript[f, p]))= (lim_(T→∞)) (1/T)|(∫_0^T)s(t) (e^(i 2π (Subscript[f, p]) t)) dt  ((|)^2)
                    ≈ (1/(M Δt))|(∫_0^(M Δt))s(t) (e^(i 2π (Subscript[f, p]) t)) dt  ((|)^2)
                   ≈ (1/(M Δt))|Δt (Sqrt[M]) FFT((Subscript[{(Subscript[s, k])}, k=1]^M) ;  p) ((|)^2)
                     ≈ ((((Δt )^2)M)/(M Δt))| FFT((Subscript[{(Subscript[s, k])}, k=1]^M) ;  p) ((|)^2)
			         =Δt| FFT((Subscript[{(Subscript[s, k])}, k=1]^M) ;  p) ((|)^2)
Now we can see why we did what we did: taking the FFT of that data and square its amplitude. Note that the quantity Δt give the unit of 1/(Whatever frequency). If t is time, it is frequency; if t is space, it is spatial frequency, we just call it generally frequency but keep in mind that it doesn't have to be strictly temporal frequency.
So the unit of the PSD is signal^2/frequency. That's why it's call power density.
Some example: if s is E is electric field, |s((|)^2) ≡|E(|^2) ≺ Intensity; if s is electric current as function of time, it is (A^2)/Hz.

Example routine

FunctionPSD[x_,DelT_]:=Module[{mx,XFFT,XPSD,fr},

(* first, we take the mean to remove the DC bias *) 
                    mx=Mean[x]  ;
(* then we take the FFT of the given series *)                    
                  XFFT=Fourier[x-mx] ;
(* PSD is simply as defined above *) 
(* but for real series, we need to take only the *)
(* positive frequency component, which is the  *)
(* first half of the FFT series *)               XPSD=Abs[Take[XFFT,Floor[Length[XFFT]/2] ] ]^2*DelT ;
(* To correspond the PDS with the frequency sries*)
(* we generate the frequency array *)                   fr=(Range[Length[XPSD]]-1.)/(DelT*Length[XFFT]) ;
(* Then, we pair the frequency with its PSD *)
(* and drop the zero frequency component *)
Drop[ Transpose[{fr,XPSD}] ,1]
                    ] ; 
LogPSD[x_,DelT_]:= Module[{fr,psd },
        {fr,psd}=Transpose[FunctionPSD[x,DelT]];
        Transpose[{fr,10*Log[10.,psd] }]
                          ] ;
PlotLogPSD[x_,DelT_]:=Module[{lpsd},
       lpsd=LogPSD[x,DelT];
       ListPlot[ lpsd
		    ,Joined->True, 
	PlotStyle->{RGBColor[0,0,1],Thickness[0.002]}
		 , PlotRange->All 
		   (* , AxesOrigin->{0,-100} *)
		, Frame->True, ImageSize->{450,300}, GridLines->Automatic] ]

LogPSDsig1=LogPSD[sig1,(1/(10^6))];
ListPlot[LogPSDsig1,Joined->True,PlotStyle->{RGBColor[0,0,1],Thickness[0.002`]},PlotRange->{-80,-20},AxesOrigin->{0,-100},Frame->True,ImageSize->{650,300},GridLines->Automatic]

We see that the signal has a component at ~ 20 kHz and at ~ 111 kHz. What is the noise power relative to the signal power? ~ 26 dB and 18 dB, respectively.
The difference between the two strong signals are ~ 8 dB. Does that make sense?
Let's see the way the signal generated:
sig1 = x1+Table[0.5*Sin[x*0.7] + 1.3 *Sin[x*0.13],{x,0,1999,1} ];
		ListPlot[sig1  , PlotJoined->True]
Strong signal amplitude = 1.3, weak signal amplitude =0.5

10*Log[10.,(1.3/0.5)^2 ]

Indeed it is about 8 dB. What about frequency? we assume the sampling interval is 10^(-6), say sec or microsecond. The frequencies are:

0.7/(2*Pi)*10^6
0.13/(2*Pi)*10^6

Which is indeed the values we observed. How about our estimation of noise?

Table[ LogPSDsig1[[i]],{i,50,60}]

Table[ LogPSDsig1[[i]],{i,800,810}]

The noise looks like the same everywhere: this is what we call white noise, because like white light, it is constant for any color except at the two signals.
We can average them over a range, say, from point 500 to 900:

rpsd=Mean[ 10.^(Transpose[
					Take[ LogPSDsig1,{500, 900}]
					                                             ] [[2]]/10.)
		                   ]

10*Log[10.,rpsd]

The noise PSD is about -60 dB/Hz. Where does this number come from?
x1=RandomArray[NormalDistribution[0,1],2000];
This expression tells us that the x1 noise is a normal (Gaussian) distribution with a standard of deviation of 1. The sampling time is 10^(-6) second, or the bandwidth is 0.5*10^6 Hz (remember Nyquist), therefore the power density is:
((  (σ^2))/(2*(0.5*(10^6)))) = (10^(-6)) /Hz or -60 dB/Hz.
Let's redo the example above:

x1=RandomReal[NormalDistribution[0,3.27],2000];
sig1 = x1+Table[0.5*Sin[x*0.7] + 1.3 *Sin[x*0.13],{x,0,1999,1} ];
		ListPlot[sig1  , Joined->True]

LogPSDsig1=LogPSD[sig1,10^-6] ;
ListPlot[ LogPSDsig1
		    ,Joined->True, 
	PlotStyle->{RGBColor[0,0,1],Thickness[0.002]}
		 , PlotRange->{-80,-20}, AxesOrigin->{0,-100}
		, Frame->True, ImageSize->{650,300}, GridLines->Automatic]

ListPlot[ FunctionPSD[sig1,10^-6]
		    ,PlotJoined->True, 
	PlotStyle->{RGBColor[0,0,1],Thickness[0.002]}
		 , PlotRange->All 
		, Frame->True, ImageSize->{650,300}, GridLines->Automatic, AspectRatio-> 300/650]

rpsd=Mean[ 10.^(Transpose[
					Take[ LogPSDsig1,{500, 900}]
					                                             ] [[2]]/10.)
		                   ]

10*Log[10.,rpsd]

compare with 
((σ^2)/(10^6)) = ((3.27^2)/(10^6))

10*Log[10., ((3.27^2)/(10^6))]

Summary

In general, it can be proven that for a Gaussian process of standard deviation σ, and a sampling interval of Δt, or a bandwidth BW=1/(2 Δt), the PSD is = (σ^2) Δt = ((σ^2)/(2 BW ))
or (σ^2) = 2 BW PSD.
The noise is taken as the square root of the PSD:
Noise amplitude = (Sqrt[Noise PSD]), and its unit is whatever/(Sqrt[Hz]).

Homework and exercise on noise and signal

Utility package - please execute this

FunctionPSD[x_,DelT_]:=Module[{mx,XFFT,XPSD,fr},

(* first, we take the mean to remove the DC bias *) 
                    mx=Mean[x]  ;
(* then we take the FFT of the given series *)                    
                  XFFT=Fourier[x-mx] ;
(* PSD is simply as defined above *) 
(* but for real series, we need to take only the *)
(* positive frequency component, which is the  *)
(* first half of the FFT series *)               XPSD=Abs[Take[XFFT,Floor[Length[XFFT]/2] ] ]^2*DelT ;
(* To correspond the PDS with the frequency sries*)
(* we generate the frequency array *)                   fr=(Range[Length[XPSD]]-1.)/(DelT*Length[XFFT]) ;
(* Then, we pair the frequency with its PSD *)
(* and drop the zero frequency component *)
Drop[ Transpose[{fr,XPSD}] ,1]
                    ] ; 
LogPSD[x_,DelT_]:= Module[{fr,psd },
        {fr,psd}=Transpose[FunctionPSD[x,DelT]];
        Transpose[{fr,10*Log[10.,psd] }]
                          ] ;
PlotLogPSD[x_,DelT_,Style_]:=Module[{lpsd},
       lpsd=LogPSD[x,DelT];
       ListPlot[ lpsd
		    ,Joined->True, 
	PlotStyle->Style
		 , PlotRange->All 
		   (* , AxesOrigin->{0,-100} *)
		, Frame->True, ImageSize->{450,300}, GridLines->Automatic] ];
		
PlotLogLogPSD[x_,DelT_,Style_]:=Module[{},
       ListLogLogPlot[ FunctionPSD[x,DelT]
		    ,Joined->True, 
	PlotStyle->Style
		 , PlotRange->All 
		   (* , AxesOrigin->{0,-100} *)
		, Frame->True, ImageSize->{450,300}, GridLines->Automatic] ]

Use your microphone input as signal source

Do the following to simulate signal measurement vs noise using your computer microphone. Play the sound below.

Sig1=SystemDialogInput["RecordSound"]

xdat1=Sig1[[1]][[1,1]] ;
DelT=1/11025  ;
PSDSig1=PlotLogPSD[xdat1,DelT,{Blue}] ;
Show[PSDSig1]

PlotLogLogPSD[xdat1,DelT,{Red}]

Now, recording some other sounds with noise

Sig2=SystemDialogInput["RecordSound"]

Perform spectral analysis: use the power spectral density function

Now you should compare the different signals - especially their noise level, signal level, and signal to noise ratio.

xdat1=Sig1[[1]][[1,1]] ;
xdat2=Sig2[[1]][[1,1]] ;
DelT=1/11025  ;
PSDSig1=PlotLogPSD[xdat1,DelT,{Blue}] ;
PSDSig2=PlotLogPSD[xdat2,DelT, {Red}];
Show[PSDSig1, PSDSig2, PlotRange-> {{1500,2500},All}]

Signal and noise

Generate some strong noise. Record it and perform spectral analysis.
Generate a signal, something you are familiar with, starting weak. Then increasing noise the signal until you see it above the noise. Determine the NES: noise equivalent signal.

2. Shot noise and related

2.1 Introduction

Suppose we measure the current from a detector or sensor or whatever, and it has a constant dark current (or DC current) even in the absence of the stimulus, (Subscript[i, dark]). We sample at regular interval Δt. The average number of electrons is:
                                    (Subscript[n, ave])  =  (Subscript[i, dark]) Δt/e
where e is the charge. The number of electrons per unit interval obeys Poisson distribution:
                                   P[n] = (((Subscript[n, ave])^n)/(n!))(e^(-(Subscript[n, ave])))
For large (Subscript[n, ave]), this distribution approx to Gaussian:
                                    P[n] ≈ (1/((Sqrt[2π]) σ))(e^(-(((n-(Subscript[n, ave])))^2)/2 (σ^(2   ))))
where σ = (Sqrt[(Subscript[n, ave])]).
What is the power spectral density of the dark current?

Exercise and Example Poisson

Select a number between 1 and 10. Use it as naverage and enter.

naverage=3.7 ;
Xpoisson=RandomInteger[PoissonDistribution[naverage] ,100]  ;
Histogram[Xpoisson]

Manipulate[ListPlot[Take[Xpoisson,n], Joined-> True],{n,1,Length[Xpoisson],1}]

Select a number between 40 and 80. Use it as "naverage" and enter.

naverage=50 ;
Xpoisson=RandomInteger[PoissonDistribution[naverage] ,2000]  ;
Histogram[Xpoisson]

Manipulate[ListPlot[Take[Xpoisson,{n,n+200}], Joined-> True, PlotRange-> {0,100}],{n,1,Length[Xpoisson]-200,1}]

Exercise and Compare Poisson to Gaussian for large mean

Compare Poisson and Gaussian distribution for large mean value: they become similar to each other. Pick a mean value >50.

mu=50.  ;
Xpoisson=RandomInteger[PoissonDistribution[mu] ,5000]  ;
XNorm=RandomReal[NormalDistribution[mu, Sqrt[mu]] ,5000]  ;  
Show[GraphicsArray[{Histogram[Xpoisson, ChartStyle-> {Red}]  
,Histogram[XNorm, ChartStyle-> Blue] }]]

Do you notice that they become very similar to each other? You can plot them on top of each other as follow:

Histogram[{Xpoisson, XNorm}]

Exercise and Dark current example 1

Consider a PD with 1 mA dark current. Suppose we measure with 20 GS/s, which means that we perform 20,000,000,000 measurements per second, what is the distribution of the electrons in each measurement like?
The mean is:                   (Subscript[n, ave])  =  (Subscript[i, dark]) Δt/e
                                         Δt=(1/(2×(10^10)))
                                        (Subscript[n, ave])  =  (10^(-3)) /(2×(10^10) e)

(10^(-3)) /(2×(10^10) ElectronCharge)

mu=312075  ;
Xpoisson=RandomInteger[PoissonDistribution[mu] ,5000]  ;Histogram[Xpoisson, ChartStyle-> {Red}]

What is the noise fluctuation?
                                σ=(Sqrt[312075])

σ=(Sqrt[312075.])

The noise current is:

σ*ElectronCharge*2×(10^10)/Sec

Hence, the current fluctuation is in the order of 1.79 microAmp. 
Is this large or small? well, it's all relative.
Consider that we have 0.1 uW signal falling on a detector with 0.5A/W responsivity. The signal would be:
                                            (Subscript[i, sig])=R P=0.05 μA
This is smaller than the noise.

Dark current exercise 2

What is the minimum optical power you can detect with 20 GS/s? (in 50 ps) with the above detector?
Repeat the same calculation as above with sampling rate at 1 MS/s. (sampling time 1 μs)
What can you conclude about noise vs. bandwidth.

Noise current exercise 3

Simulate the current (including noise of course) for a detector with 5 μA dark current for the following cases:
1- sampling at 10 GS/s (10 gigasamples per second)
2- sampling at 100 MS/s (100 megasamples per second)
In each case, calculate the number of average electrons you detect, their fluctuation (standard of deviation). Which case has larger electron (hence current) fluctuation? If so, which case (1 or 2) is noisier? (BE VERY THOUGTHFUL and CAREFUL about this question).

2.2 Shot noise

Even without dark current, any signal also has statistical fluctuation that obeys Poisson statistics:
                                                (Subscript[n, ave])  =  (Subscript[i, signal]) Δt/e
                                               P[n] = (((Subscript[n, ave])^n)/(n!))(e^(-(Subscript[n, ave]))) 
The corresponding distribution of the current measurement is:
                                   P[i] ≈ (1/((Sqrt[2π]) (Subscript[σ, i])))(e^(-(((i-(Subscript[i, signal])))^2)/2 ((Subscript[σ, i])^(2   ))))
where (Subscript[σ, i]) =(Sqrt[(Subscript[n, ave])])(e/(Δt)). Notice that:
                                       ((Subscript[σ, i])^2) =(((Subscript[n, ave]) (e^2))/((Δt)^2))=  (Subscript[i, signal]) (e/(Δt)) =(Subscript[i, signal]) e 2B
Or:                                            (OverBar[(Subscript[i, noise]^2)])=2 e (Subscript[i, signal]) B  
(remember that 1/Δt = 2B Nyquist); then PSD of this distribution (number of electrons in Δt interval) is:
                                          (((Subscript[σ, i])^(2 ))/(2 B))= (((Subscript[i, signal]) e 2B)/(2 B)) =  (Subscript[i, signal]) e 
Very often, for photodetector, the dark current shot noise is the limiting factor, then people define a noise equivalent power, which is the power that generates a signal = noise.
                                            NEP = ((Sqrt[2  (Subscript[i, dark]) e B])/R)
The unit of NEP is __W__
Often, we also scale the NEP per unit of (Sqrt[B]), hence the NEP/(Sqrt[B]) is:
                                            NEP = ((Sqrt[2  (Subscript[i, dark]) e])/R)
The unit of NEP per bandwidth is __W/(Hz^(-1/2))_____
and the detector detectivity is define as:
                                            D = (1/NEP) = (R/(Sqrt[2  (Subscript[i, dark]) e]))
For semiconductor PD, the dark current is proportional to the detector area, thus, people define a specific detectivity by scaling D vs. detector area (Subscript[A, d]):
                                          (SuperStar[D]) = (1/(SuperStar[NEP])) = (R/(Sqrt[2 (Subscript[i, dark]) e/(Subscript[A, d])]))= R (Sqrt[((Subscript[A, d])/(2  (Subscript[i, dark]) e))]) 
The unit of specific detectivity (SuperStar[D])  is thus:   cm (Hz^(1/2))/W.
Notice, in some case, you may see this (without the factor of 2). It is just a matter of convention about the noise power (just one band or both bands).
                                    (SuperStar[D]) = (1/(SuperStar[NEP])) = (R/(Sqrt[ (Subscript[i, dark]) e/(Subscript[A, d])]))= R (Sqrt[((Subscript[A, d])/( (Subscript[i, dark]) e))])

Use the formulas given above to calculate the detectivity for a detector with a dark current of 0.5 nA, and responsivity of 0.65 A/W

If you are trying to detect a 10-nW signal, and suddenly, someone turns the light on in the room with 10 mW of ambient light falling on the detector, what is your SNR?

2.3 Optical quantum noise

Similar to electronic shot noise, we can have statistical fluctuation with photons. However, photon statistics is more complex when they are in certain quantum states. Neglecting these cases of quantum coherence, we can treat most other cases with Poisson statistics.
Suppose we have a stream of photons with an expected average (Subscript[n, ave]) within some time interval Δt. Then probability of detecting n photons is:
                                                    P[n] = (((Subscript[n, ave])^n)/(n!))(e^(-(Subscript[n, ave])))
But this is what we detect, it is NOT necessarily the statistics of light from a source, since we don't know the source characteristics, i. e. how it emits light.
For quantum-incoherent light emitters, the intensity distribution is actually:
                                                    P[I] = (1/(Subscript[I, ave]))(e^(-I/(Subscript[I, ave])))
Thus, the photon statistics from such a source is obtained by averaging over the Poisson distribution, which yields:
                                                  P[n] = (((Subscript[n, ave])^n)/(((1+(Subscript[n, ave])))^(n+1)))     
This is known as Bose-Einstein distribution.

Example, let nave=9, this is the distribution:

BarChart[Table[(((Subscript[n, ave])^n)/(((1+(Subscript[n, ave])))^(n+1))) /. (Subscript[n, ave])-> 9, {n,0,50}]]

Here is with larger nave.

BarChart[Table[(((Subscript[n, ave])^n)/(((1+(Subscript[n, ave])))^(n+1))) /. (Subscript[n, ave])-> 100, {n,0,500}] ]

2.4 Theoretical digital signaling quantum limit

Discussion

If we have a perfect detector, one that gives us zero dark current, and 100% quantum efficiency, and we use it for optical communication, do we have the probability of errors?
The answer is yes! this has NOTHING to do with how perfect the detector is, but the fundamental quantum process of detection. Remember from above, suppose we have a pulse of light with an expected average (Subscript[n, ave]), we don't actually detect (Subscript[n, ave]), but can be any number n with the probability:
                                                       P[n] = (((Subscript[n, ave])^n)/(n!))(e^(-(Subscript[n, ave])))

BarChart[Table[(((Subscript[n, ave])^n)/(n!))(E^(-(Subscript[n, ave]))) /. (Subscript[n, ave])-> 4, {n,0,12}] ]

We see that there is a finite probability of error: we may detect no photon while in fact there is an average (Subscript[n, ave])>0:
                                              P[0] = (((Subscript[n, ave])^0)/(0!))(e^(-(Subscript[n, ave])))=(e^(-(Subscript[n, ave]))) 
Suppose our system is digital, so that if there is one photo-electron, we call bit 1 and without it, we call bit 0, then in fact we have an error when we detect 0 photoelectrons. In other words, we miss a bit 1 signal.
Hence the probability of error is:
                                            P[e|bit 1] = (((Subscript[n, ave])^0)/(0!))(e^(-(Subscript[n, ave])))=(e^(-(Subscript[n, ave])))

Bit error rate

We define bit-error-rate (BER) is the fraction of error of all the bits. Hence, a BER of (10^(-9)) means that out of one billion bits, we may expect 1 bit error. Notice that we say "expect" only, not guarantee that there must be one bit error. We may have more or less, but the ensemble average is 1 bit of error in one billion bits.
Now, given the error above, what is the BER as a function of (Subscript[n, ave])?

Plot[Log[10,(E^(-(Subscript[n, ave]))) ],{(Subscript[n, ave]),1,30}, Frame-> True, GridLines-> Automatic, AspectRatio-> 1.5]

So, for a BER of (10^(-9)), we need at least:
                               (Subscript[n, ave])=-Log[BER]

-Log[10.^-9]

Or 20.7 detected photons. Of course, this is the quantum limit for ideal receiver. Realistic system is a lot worse than that as we will see.

Example: For 10 Gb/s and a wavelength 1.55 um, how much power we need for an error rate of 10^-9?

The number of photons per bit:
                                  (Subscript[n, ave])=-Log[BER]
The number of photons per second:
                                   (Subscript[N, p])=(Subscript[n, ave])BR=-Log[BER] BR
 The power is:              P=(Subscript[N, p]) h ν=(Subscript[n, ave])BR h ν=-Log[BER] h ν BR

<<PhysicalConstants`;
λ=1.55 ; BER=10^-9 ; BR=10^10 ;
-Log[BER] (1.2398419/λ)  ElectronCharge *Watt/Coulomb BR

A mere 27 nWatt. In dBm, it is:

10*Log[10, 2.6558476768660435`*^-8 *10^3]

A useful physical constant is:

PlanckConstant *SpeedOfLight

Or -46 dBm.
A chart we often plot is the BER as a function of power for a given BR:
                                     BER=  (P/((h ν) BR)) =((P λ)/(h c BR))

BERf[P_,λ_, BR_] :=Exp[-((P  λ)/(1.9864455003959036`*^-19  BR)) ]

BERf[10^(-45.76/10)*10^-3,1.55,10^10]

As we can see, we obtain the value of (10^(-9)) BER as expected.

Plot[Table[Log[10., BERf[10^(P/10)*10^-3,1.55,10^br]],{br,8,10,1}]
, {P,-70,-40}, PlotRange-> {-15,-3}, PlotStyle-> {Red, Green, Blue}, Frame-> True ,GridLines-> Automatic, AspectRatio-> 1.5 ]

2.5 Theoretical analog signaling quantum limit

Comparison of important parameters between DIGITAL and ANALOG
                             DIGITAL                              ANALOG
                              BER                                     Signal-to-noise ratio (SNR, S/N)
                              photons/bit                           Photocurrent
                              Discrete signal (1,0)             Continuous
                              Bit rate                                 Bandwidth (analog)                                  
The signal of analog system is not bit, but signal current (continuous), (Subscript[i, signal]). The noise is                             
                                                         (OverBar[(Subscript[i, noise]^2)])=2 e (Subscript[i, signal]) B 
The signal to noise ratio (SNR) in terms of current is:
                                                    (Subscript[SNR, curr])=((Subscript[i, signal])/(Sqrt[2 e (Subscript[i, signal]) B ]))=(Sqrt[((Subscript[i, signal])/(2 e B ))]) 
There is another definition of SNR that is in terms of signal power, not current.
Remember that Power ≺ (Current^2):
                                                  (Subscript[SNR, pow])=((Subscript[i, signal])/(2 e B )) 
In terms of incident optical power:
                                               (Subscript[SNR, pow])=((Subscript[i, signal])/(2 e B )) =η(P/(h ν)) e (1/(2 e B)) =η(P/(2 h ν B))     
Obviously, we are assuming an ideal detector without dark current. Assuming also that the QE is perfect 1,
                                              (Subscript[SNR, pow])=(P/(2 h ν B))

Example plot of signal to noise ratio

SNR[P_,λ_,f_]:= ((P  λ)/(2*f*1.9864455003959036`*^-19 ))  ;

Plot[Table[10*Log[10,SNR[10^(P/10-3),1.55,10^bw] ],{bw,6,10,1}]
,{P ,-60,-20}
,Frame-> True, GridLines-> Automatic, AspectRatio-> 1.5]

3. Other noise

3.1 Johnson (thermal) noise

Electrons in a medium at a finite temperature do not rest, but agitate. Thus, there is a fluctuating current, of which the average is zero, but not the mean square. The kinetic energy of electrons ~ ∝(Subscript[k, B])T. The noise power is:            (Subscript[P, thermal])=2(((Subscript[k, B])T)/(Δt))=4 (Subscript[k, B])T B
If the curcuit resistor is R, then             (i^2)R =4 (Subscript[k, B])T B      
Or:                                                           (i^2) =((4 (Subscript[k, B])T B)/R)

<<PhysicalConstants`;
BoltzmannConstant

(Sqrt[((4 (Subscript[k, B])T B)/R)])  /. (Subscript[k, B])-> BoltzmannConstant /. T-> 300 Kelvin  /. R-> 50 Ohm /. B-> 10^6 Hz

The noise current is ~ 18 nA.

3.2 Relative intensity noise and laser noise

3.2.1 Introduction

Beyond the light source quantum noise, lasers have other noise sources that come from a number of complex physical mechanisms. One way to describe the laser noise is the concept of relative intensity noise or RIN. As you measure the laser power, for example 1 mW, you will see that it is NOT always exactly 1 mW but fluctuate a little bit. See below.

xrin=RandomReal[NormalDistribution[100,0.3],1000] ;
ListPlot[xrin, Joined-> True, PlotRange-> {0,105}]

ListPlot[xrin, Joined-> True, PlotRange-> {98,102}]

If we take the fluctuation, divide it by the average, we will see the change in percentage. The above curve shows a fluctuation ~ 1%.

3.2.2 General discussion

General concept: Suppose we have a signal S[t] that is supposed to be constant (Subscript[S, 0]). We measure the fluctuation in terms of its mean value (Subscript[S, 0]): n[t]=((S[t]-(Subscript[S, 0]))/(Subscript[S, 0])). We call this relative amplitude noise. If another quantity is a power function of S[t], e. g. P[t]=((S[t])^r), (Subscript[N, p])[t]=((P[t]-(Subscript[P, 0]))/(Subscript[P, 0]))= ((((S[t])^r)-((Subscript[S, 0])^r))/((Subscript[S, 0])^r))≈(((S[t]-(Subscript[S, 0]))r ((Subscript[S, 0])^(r-1)))/((Subscript[S, 0])^r))≈ r n[t]. Thus the power spectral density of relative noise of P[t] is (r^2) of that of S[t].

Back to laser RIN: a laser source intensity (or power) is not constant, but fluctuates. Sure, we see that at the very least, we have quantum noise fluctuation. But usually, the fluctuation is even much larger than quantum noise. They are called intensity noise. Since it is measured relative to the DC intensity level, it is called RIN: Relative intensity noise
We express the power as: P(t), and let (Subscript[P, 0]) = <P(t)>. The relative intensity signal is: r(t) = ((P(t))/(Subscript[P, 0]))(we can drop 1, which is just a constant with zero noise). Its noise PSD is:
       RIN(f)= (lim_(T→∞)) (1/T)|(∫_0^T)r(t) (e^(i 2π f t)) dt  ((|)^2)
Since r(t) is unitless, the PSD unit is usually expressed as dB/Hz. Sometimes, we just use the square root value and still call the same name RIN. How to tell one from the next? look at the unit, whether dB/Hz (power) or dB/(Sqrt[Hz]) (amplitude).

Understanding RIN is very important to use the laser properly.

3.2.3 Typical laser RIN

Here is ONLY an example. By no means the same for all lasers. Each laser has its own RIN PSD.

Example: Suppose you purchase a laser diode with a flat RIN of -125 dB/Hz from 10-100 MHz. You want to modulate the laser at 50 MHz. What is the intensity noise do you expect?
                                               (Noise^2)=RINPSD × BW  
                                               (Noise^2)=(10^(-12.5))/Hz × 50 MHz =(10^(-12.5))/Hz × 5×(10^7) Hz 
                                               (Noise^2)=5 ×(10^(-5.5))
                                              Noise=(Sqrt[5 ×(10^(-5.5))])

N[(Sqrt[5 *(10^(-5.5))]) ]

Thus, the relative noise is 0.3976% or ~ 0.4%. This is how much the laser fluctuates at 50 MHz.
If you modulate the laser at 12.5 MHZ instead, what is the fluctuation?

4. Receiver (6323 only - part 2)

Introduction

Link to ppt file - receiver

Response

With voltage amplifier the response is simply the net load resistance, which is parallel between the bias load and the amplifier load
                                   (Subscript[R, net load])=(((Subscript[R, bias]) (Subscript[R, amp]))/((Subscript[R, bias]) +(Subscript[R, amp])))

With transimpedance amplifier:
                                  H[f] ~  ((  -(Subscript[R, f]))/(1+ i 2 π f (Subscript[R, f]) (Subscript[C, T]) /G  )) 
where G is the open-loop gain of the op-amp.
The approx 3-dB bandwidth is:
                                        2 π f (Subscript[R, f]) (Subscript[C, T]) /G ~ 1
                                        or: f~ (G/(2 π  (Subscript[R, f]) (Subscript[C, T])))
Notice that gain (given by (Subscript[R, f])  ) is trade-off with bandwidth. This gain-bandwidth trade-off is very fundamental. Recall the gain-BW product of transistor. Later on we will also see gain-BW product limit of APD.

For other circuits, e. g. p-i-n FET (HEMT) or p-i-n HBT or APD-FET, HBT etc... all can also be modelled generically   H[f] ~  ((P[f])/(Q[f]))  where P and Q are polynomials. Often, parameters are obtained empirically.

Noise figure (general discussion)

A very important concept associated with any systems is noise figure. It is most often associated with amplifiers. Suppose we have a weak signal and we want to amplify.
Suppose we have this signal:

x1=RandomReal[NormalDistribution[0,1],2000];
sig1=(1/(10^2))(0.1` x1+Table[0.5` Sin[x 0.7`]+1.3` Sin[x 0.13`],{x,0,1999,1}]);
s1p=ListPlot[sig1,Joined->True,PlotStyle->{RGBColor[1.`,0.`,0]}]

LogPSDsig1=LogPSD[sig1,(1/(10^6))];
psd1p=ListPlot[LogPSDsig1,Joined->True,PlotStyle->{RGBColor[1,0,0],Thickness[0.002`]},PlotRange->{-140,-60},AxesOrigin->{0,-100},Frame->True,ImageSize->{650,300},GridLines->Automatic]

Our amplifier has a power gain of 30 dB. After amplification, we have a stronger signal:

x2=RandomReal[NormalDistribution[0,1],2000];
sig2=(10^(1.5`)) sig1+0.1` x2;
s2p=ListPlot[sig2,Joined->True,PlotStyle->{RGBColor[0.`,0.`,1]}]

We see that our signal is stronger. It is ~ 0.45 V compared with 0.015 V before, a gain of ~ 15 dB (power gain is 30 dB). Let's compare the signals. Do you think the noise is less for the amplified signal?

Show[s1p, s2p]

We can blow up the first signal to compare with the amplified signal

s1large=ListPlot[sig1 30,Joined->True,PlotStyle->{RGBColor[1.`,0.`,0]}]

Show[ s2p,s1large]

Why aren't the signals exactly on top of each other?
Let's look at the amplifed signal power

LogPSDsig2=LogPSD[sig2,(1/(10^6))];
psd2p=ListPlot[LogPSDsig2,Joined->True,PlotStyle->{RGBColor[0,0,1],Thickness[0.002`]},PlotRange->{-100,-30},AxesOrigin->{0,-100},Frame->True,ImageSize->{650,300},GridLines->Automatic]

It looks quite similar to the one before amplified, with the only exception that the signal power is higher, ~ -43 dB instead of -73 dB.
Let's compare them on the same scale

Show[psd1p, psd2p, PlotRange->{-140,-40}, GridLines->Automatic]

Are they similar with each other except for a vertical shift?

We notice that the signal gain 30 dB in power, but the noise floor seems to gain 40 dB in power. Thus, the signal-to-noise ratios before, 47 and 40 dB, become only ~ 37 and 30 dB. What happens? the SNR is worse! This happens because no amplifier is perfect. All amplifers add some extra noise to the output signal. The degradation of SNR is called noise figure (NF). Noise figure thus has a very specific definition, it's the noise level added to a signal by the amplifer. It is a figure of merit of the amplifier; needless to say, the lower it is the better.

If amplifers add noise, why bother using them?

Noise model p-i-n

Additive noise model

We will follow the additive noise model, which states that the total noise is the sum of the square of uncorrelated noise, which is based on the Gaussian:
                                             (Subscript[i, total]^2)=(Subscript[i, n1]^2)+(Subscript[i, n2]^2)  + ...
For a PIN receiver, there are 3 effective noise terms: shot noise, thermal noise, and amplifier noise.
Shot noise:                                 (OverBar[(Subscript[i, sig]^2)])=2 e  ((Subscript[I, photo])+(Subscript[I, dark])) B
Thermal:                                     (OverBar[(Subscript[i, therm]^2)]) =((4 (Subscript[k, B]) T B)/(Subscript[R, L]))
For the amplifier, a most general model is that it is a source of noise that has a series voltage noise (Subscript[v, a])[f] and a shunt current noise ((Subscript[i, a])^2)[f]. Let Y be the shunt admittance of the amplifier, The total noise is an integration of these two over all frequencies within the bandwidth:
                                               (OverBar[(Subscript[i, Amp]^2)]) =(∫_o^B)(Subscript[v, a])[f] ((Subscript[i, a])^2)[f] (((|Y|))^2)df     
which is usually done with an empirical model.
The total noise is:
                               (Subscript[i, total]^2)=2 e M ((Subscript[I, photo])+(Subscript[I, dark])) B+((4 (Subscript[k, B]) T B)/(Subscript[R, L]))  + (OverBar[(Subscript[i, Amp]^2)])   
Very often the amplifer noise  (OverBar[(Subscript[i, Amp]^2)])  is referred to the input load resistance  (Subscript[R, L]) for an effective noise power of the amplifier. Then we write:
                             (Subscript[i, total]^2)=2 e M ((Subscript[I, photo])+(Subscript[I, dark])) B+((4 (Subscript[k, B]) T B)/(Subscript[R, L]))  (1+(1/(4 (Subscript[k, B]) T B))(Subscript[R, L])   (OverBar[(Subscript[i, Amp]^2)])) 
Then, the term    (1+(1/(4 (Subscript[k, B]) T B))(Subscript[R, L])   (OverBar[(Subscript[i, Amp]^2)])) is treated as a parameter that is characteristic of the amplifier which we call noise figure:
                                   (1+(1/(4 (Subscript[k, B]) T B))(Subscript[R, L])   (OverBar[(Subscript[i, Amp]^2)]))≡ (Subscript[F, n]) 
This is indeed a very common and key concept as discussed above.

SNR

We now can obtain the SNR:
                              SNR= (( (Subscript[I, photo]^2))/(B (2 e M ((Subscript[I, photo])+(Subscript[I, dark]))   + 4 (Subscript[k, B]) T (Subscript[F, n]) /(Subscript[R, L])) ))

Noise model APD

Excess noise

A most important concept with APD noise is excess noise. When gain is applied, the ideal condition is that:
                                           (Subscript[I, out])=M ((Subscript[I, photo])+(Subscript[I, dark]))
where M is the gain. In reality, the noise is higher than the shot noise expected for  (Subscript[I, out]), which is:
                                    (OverBar[(Subscript[i, ideal]^2)])=2 e M ((Subscript[I, photo])+(Subscript[I, dark])) B
We will see that in reality, no amplification or gain process is perfect: there will always be extra noise introduced (see noise figure concept later) on top of the input noise. For APD, there is intrisic random gain process (electron ionization process is a stochastic process with certain probability distribution). This process is intrinsic to the semiconductor and cannot be controlled.
Thus, a theoretical model can show that the actual noise is:
                                    (OverBar[(Subscript[i, APD n]^2)])=2 e (M^(2+x)) ((Subscript[I, photo])+(Subscript[I, dark])) B    
where x is a parameter that is intrinsic to the carrier and the semiconductor (a complex dependence on electron mass, hole mass, ionization probability... ). Typical value of x is ~ 0.2 - 1 for various semiconductors and carriers.
The SNR of APD is thus:
                                           SNR= (((M^2) (Subscript[I, photo]^2))/(B (2 e (M^(2+x)) ((Subscript[I, photo])+(Subscript[I, dark]))  + 4 (Subscript[k, B]) T (Subscript[F, n]) /(Subscript[R, L])) ))
where (Subscript[F, n])  is the noise figure.

SNRAPD[Ip_,Id_, T_, R_, Fn_,B_,M_, x_ ]:= Module[{q, kB},
q=1.6*10^(-19)  ; kB=1.3806504`*^-23   ;
M^2 Ip^2 / (B*(2 q *(Ip+Id) M^(2+x) + 4 kB * T *Fn/R))  ]

LogLinearPlot[Table[10*Log[10, SNRAPD[10^-Is,5*10^-10, 300, 10^6, 2,10^6,M,0.3 ]],{Is, 4,10,1}], {M,1,100}, AspectRatio-> 1, Frame-> True, GridLines-> Automatic]

We notice that depending on the signal, there is an optimal SNR for gain. In other words, when the signal is small, we need gain, but not as large as possible. Only certain value of M will give us the best SNR. In fact, we see that:

FullSimplify[(∂_M)(((M^2) (Subscript[I, photo]^2))/(B (2 e (M^(2+x)) ((Subscript[I, photo])+(Subscript[I, dark]))  + 4 (Subscript[k, B]) T (Subscript[F, n]) /(Subscript[R, L])) ))]

What if (-4 T (Subscript[k, B]) (Subscript[F, n])+e x ((Subscript[i, dark])+(Subscript[i, photo])) (Subscript[R, L]) (M^(x+2)))==0 ? Is it possible? Yes.

Solve[(-4 T (Subscript[k, B]) (Subscript[F, n])+e x ((Subscript[I, dark])+(Subscript[I, photo])) (Subscript[R, L]) (M^(x+2)))==0 , {M}]

FullSimplify[{{M->(4^(-(1/(-x-2)))) ((((((Subscript[I, dark]) e x (Subscript[R, L])+e (Subscript[I, photo]) x (Subscript[R, L]))/(T (Subscript[k, B]) (Subscript[F, n])))))^(1/(-x-2)))}}]

Indeed, the reason is because of the excess noise factor. Let x=0

LogLinearPlot[Table[10*Log[10, SNRAPD[10^-Is,5*10^-10, 300, 10^6, 2,10^6,M,0 ]],{Is, 4,10,1}], {M,1,100}, AspectRatio-> 2, Frame-> True, GridLines-> Automatic]

This shows that any gain is good and the larger the better!! Because with gain, we defeat the thermal noise and hence, the larger the better. But because of the excess noise factor, too high gain is counterproductive because the excess noise also grow and cause worse SNR. Infact, compare small x and large x:

ap=LogLinearPlot[Table[10*Log[10, SNRAPD[10^-Is,5*10^-10, 300, 10^6, 2,10^6,M,0 ]],{Is, 4,10,1}]
, {M,1,100}, AspectRatio-> 2, Frame-> True, GridLines-> Automatic, PlotStyle-> {Thick, Red}] ;
bp=LogLinearPlot[Table[10*Log[10, SNRAPD[10^-Is,5*10^-10, 300, 10^6, 2,10^6,M,1 ]],{Is, 4,10,1}]
, {M,1,100}, AspectRatio-> 2, Frame-> True, GridLines-> Automatic,  PlotStyle-> {Thick, Blue}] ;

Show[ap, bp]

More accurate model of excess avalanche noise factor

The excess noise factor is actually more complex than just just (M^x). It is different for electrons and holes, because they have different transport behavior and ionization rate.
A more accurate model for excess noise factor has the factor different for electrons and holes. For electron:
                                 F[M]= M(1-(1-k)(((((M -1)/M)))^2) )  ;
for hole: 
                                  F[M]= M(1+(((1-k)/k))(((((M -1)/M)))^2) ) 
where k is the ratio of ionization coefficients for electron and hole.
Obviously, k is a very important factor and one can plot

klist={0.01,0.02,0.05, 0.1,0.2,0.5,1,2,5,10,20,50,100}  ; 
LogLogPlot[Table[ M (1-(1-klist[[i]])(((((M -1)/M)))^2) ) ,{i,1, Length[klist]}], {M,1,1000}, AspectRatio-> 1, PlotStyle-> {Table[Hue[(j-1)/Length[klist],1,1],{j,1, Length[klist]}] }  ]
LogLogPlot[Table[ M (1+((1-klist[[i]])/klist[[i]])(((((M -1)/M)))^2) ) ,{i,1, Length[klist]}], {M,1,1000}, AspectRatio-> 1, PlotStyle-> {Table[Hue[(j-1)/Length[klist],1,1],{j,1, Length[klist]}] }  ]

Noise and BER with Gaussian model

If we have a string of 0 and 1 bit, noise can cause error when bit 1 is identified as 0 and vice versa. A threshold is usually set such that if signal S>(Subscript[s, T]), it's identified as 1 and 0 otherwise.
For a Gaussian noise model:

Manipulate[
If[TrueQ[RUNnoisebitdisc],
Manipulate[
lmax=Max[(1/(σ0*(Sqrt[2 π]))),(1/(σ1*(Sqrt[2 π])))];
x0=((σ0*xHi)/(σ0+σ1)) ; 
BER0=Log[10.,(1/2) (1.+Erf[((x0-xHi)/((Sqrt[2])σ1))])+(1/2) Erfc[(x0/((Sqrt[2]) σ0))]]  ;
BER=Log[10.,(1/2) (1.+Erf[((xT-xHi)/((Sqrt[2])σ1))])+(1/2) Erfc[(xT/((Sqrt[2]) σ0))]]  ;
{Show[Plot[{Exp[-x^2/(2*σ0^2)]/(σ0*(Sqrt[2 π]))
,Exp[-(x-xHi)^2/(2*σ1^2)]/(σ1*(Sqrt[2 π]))}
,{x,-0.5,2}
,PlotStyle-> {{Thickness[0.005],Red},{Thickness[0.005],Blue}}
,Filling-> {{1-> {0,Hue[0.,0.5,1.,0.7]}}
,{2-> {0,Hue[0.65,0.5,1.,0.7]}}}
,PlotRange-> {All,{0.,pmax}}
,AxesLabel-> {"Signal level","Probability density"}
,LabelStyle-> {14,FontFamily-> "Arial"}]
,Graphics[{Hue[0.36,1.,0.6],Thickness[0.007],Line[{{xT,0.},{xT,lmax}}]}]
,AspectRatio-> 0.3, ImageSize-> 500
,Epilog->{Inset[ Mouseover["Discriminant\nthreshold"
,"Move threshold for\nminimum BER (blue arrow)"],{xT,Min[0.9*pmax,1.15*lmax]}
,BaseStyle-> {10,White,Bold,FontFamily-> "Arial"}
,Background-> Hue[0.4,1.,0.3]]}]

,Graphics[{Red,Rectangle[{1.,Min[Max[BER,-15],-3]},{0,-15}],
Hue[0.35,1,0.8],Rectangle[{0.,-3},{1.,Max[Min[BER,-3],-15]}]
,White,Table[Line[{{0,m},{1,m}}],{m,-15,-3}]
,Black,{Thickness[0.025],Line[{{-0.2,-9},{1.,-9}}]}
,{Blue,Thickness[0.05],Arrowheads[0.3],Arrow[{{-2,Min[Max[BER0,-15],-3]}
,{0.,Min[Max[BER0,-15],-3]}}]}
,Table[Inset[NumberForm[m,{2,0}],{1.8,m}
,BaseStyle-> {10,Bold,FontFamily-> "Arial"}],{m,-15,-3}]
,Inset["BER (10^)",{0.25,-2.5}
,BaseStyle-> {14,Bold,FontFamily-> "Arial"}]}]}
,Grid[{{Control[{{xHi,1.},0.5,2.}]
,Control[{{σ1,0.1,"bit 1 noise"},0.05,1}]},{Control[{{pmax,5.},2.,20}],Control[{{σ0,0.1,"bit 0 noise"},0.05,1}]}
,{,Control[{{xT,0.8,"Discriminant\nthreshold"},-0.5,1.5}]
}}]
]
,Button[,RUNnoisebitdisc=True ]]
,Column[{Show[,ImageSize-> {400,35}] 
,Button["ECE 6323 - APP - 2-LEVEL DIGITAL DISCRIMINATION DETECTION",
BaseStyle->{18,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]}]
,Row[{Button["RUN STATUS/CONTROL",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUNnoisebitdisc,False,"->"},{False,True},ControlType-> Checkbox}] 
,Spacer[35]
,Button[Dynamic[If[TrueQ[RUNnoisebitdisc],,]]
,{ If[TrueQ[RUNnoisebitdisc],RUNnoisebitdisc=False,RUNnoisebitdisc=True]}]}]
]

Pbit0[x_,noise_]:=Exp[-x^2/(2*noise^2)]/(noise*(Sqrt[2 π]))  ;
Pbit1[x_,noise_]:=Exp[-(x-1)^2/(2*noise^2)]/(noise*(Sqrt[2 π]))  ;
Plot[{Pbit0[x,0.1],Pbit1[x,0.1]},{x,-0.4,1.4}, PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}}
, PlotRange-> All, Filling-> Axis ]

We see that it seems easy to distinguish bit 1 from bit 0 here.

Pbit0[x_,noise_]:=Exp[-x^2/(2*noise^2)]/(noise*(Sqrt[2 π]))  ;
Pbit1[x_,noise_]:=Exp[-(x-1)^2/(2*noise^2)]/(noise*(Sqrt[2 π]))  ;
Plot[{Pbit0[x,0.25],Pbit1[x,0.25]},{x,-0.5,1.5}, PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}}
, PlotRange-> All , Filling-> Axis]

Pbit0[x_,noise_]:=Exp[-x^2/(2*noise^2)]/(noise*(Sqrt[2 π]))  ;
Pbit1[x_,noise_]:=Exp[-(x-1)^2/(2*noise^2)]/(noise*(Sqrt[2 π]))  ;
Plot[{Pbit0[x,0.25],Pbit1[x,0.1]},{x,-0.5,1.5}, PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}}
, PlotRange-> All , Filling-> Axis]

We see now that we really have a problem if the signal occurs near center. We can choose the threshold at 0.5. Everything above is 1 and 0 otherwise. But there is ambiguity. We see that there is a significant probability that we can be wrong. That's the tail end of both distribution.

We can calculate the error rate. This is the error for mistaking 0 as 1 (because the noise make it > 0.5)

(1/(σ(Sqrt[2 π])))(∫_(1/2)^∞)(E^(-((x^2)/(2  (σ^2)))))\[DifferentialD]x

Or: (Subscript[P, e])[1|0]=(1/2) erfc((1/( 2 σ (Sqrt[2]))))  
This is the error for mistaking 1 as 0 (because the noise make it < 0.5)

(1/(σ(Sqrt[2 π])))(∫_(-∞)^(1/2))(E^(-((((x-1))^2)/(2  (σ^2)))))\[DifferentialD]x

Or: (Subscript[P, e])[0|1]=(1/2) erfc((1/( 2 σ (Sqrt[2]))))   
 Not surprisingly, they are the same because of the symmetry across the line 0.5.
 Total error rate is: 
                                                   BER=erfc((1/( 2 σ (Sqrt[2]))))

BERGauss[σ_]:=erfc((1/( 2 σ (Sqrt[2]))))  ;
LogPlot[BERGauss[σ], {σ,0.07,0.2}, AspectRatio-> 1.5, Frame-> True, GridLines-> Automatic]

We can express in terms of signal to noise ratio:

BERGauss2[SNR_]:=erfc((SNR/( 2  (Sqrt[2]))))  ;
LogPlot[BERGauss2[10^(snr/10.)], {snr,1,12}, AspectRatio-> 1.5, Frame-> True, GridLines-> Automatic]

Often, we plot the BER in the reverse log scale:

BERGauss2[SNR_]:=erfc((SNR/( 2  (Sqrt[2]))))  ;
LogPlot[-Log[10,BERGauss2[10^(snr/10.)]], {snr,2,11}, AspectRatio-> 1.5, Frame-> True, GridLines-> Automatic, PlotRange-> {1,12}]

BERGauss2[SNR_]:=erfc((SNR/( 2  (Sqrt[2]))))  ;
Plot[-Log[-Log[10,BERGauss2[10^(snr/10.)]]], {snr,2,11}, AspectRatio-> 1.5, Frame-> True, GridLines-> Automatic]