Numerical Methods for Fourier Analysis 
of LTI
- part 2

ECE 3336/3337/3340/4375

Han Q. Le (c)        Also for ECE6323

6. Power spectral density

6.1 Intro (for review if no prior PSD discussion )

Consider a phasor like the App: what is the power across a load resistor? We learn that
                                      P[t] = v[t] i[t] = (R i[t]) i[t]= R ((i[t])^2)  
In phasor term, we take the real part of a physical quantity, hence:
                                       P[t] = R ((Re[i[t]])^2) = R ((Subscript[i, 0])^2) (((Re[(e^(i (ω t+ϕ)))]))^2) =R ((Subscript[i, 0])^2) (Cos^2)[ω t+ϕ]

Clearly, we are not really interested in the harmonic time-variation of power. The only thing matters is the time-average power:
                                  ⟨P⟩ =(1/T)(∫_0^T)P[t]dt=R ((Subscript[i, 0])^2) (1/T)(∫_0^T)(Cos^2)[ω t+ϕ]dt 
 By choosing T to be a multiple of the period:
                                                 ⟨P⟩ =(1/2)R ((Subscript[i, 0])^2) 
  As it turns out, this is a very generic result for all harmonic quantity s[t]:
                                 (((Re[s[t]]))^2)=(1/4)(s[t]+(SuperStar[s])[t]) (s[t]+(SuperStar[s])[t])  
                                                  =(1/4)(((s[t])^2)+(((SuperStar[s])[t])^2)+ 2 s[t](SuperStar[s])[t] )

Since the time-average of ((s[t])^2) and (((SuperStar[s])[t])^2) are zero, because of harmonic behavior, it becomes:
                                              ⟨((OverBar[s])^2)⟩ =(1/2)((|s|)^2)

What happens if the signal consists of many harmonic signals?
                                              s[t]= (∑_(m=1)^∞)(Subscript[a, m]) (e^(i ((Subscript[ω, m])t+(Subscript[ϕ, m]))))

As it turns out, this principle is generally true because all the cross harmonic terms also vanish over infinite time:
                                             ⟨((OverBar[s])^2)⟩ =(Underscript[Lim, T-> ∞])(1/T)(∫_0^T)((| (∑_(m=1)^∞)(Subscript[a, m]) (e^(i ((Subscript[ω, m])t+(Subscript[ϕ, m]))))|)^2)  dt 
                                                     =(1/2) (∑_(m=1)^∞)((|(Subscript[a, m])|)^2)

6.2 Power spectral density (for review if no prior PSD discussion)

The meaning of each (1/2)((|(Subscript[a, m])|)^2) is that it is the power of the spectral component of frequency (Subscript[ω, m]) of the signal. The signal total power is a sum of those of all spectral components.

What if s(t) is a general empirical signal that we don’t know its analytic form?
Then we define the power-spectral-density (PSD) function of a signal s(t) as:  
                            PSD(f)= (Underscript[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 
       PSD(f)=(Underscript[Lim, T→∞])  (1/T)|(∫_0^T)s(t) (e^(i 2π f t)) dt ((|)^2)= (Underscript[Lim, T→∞]) (1/T)|(∫_0^T)A ((e )^(-i 2 π (Subscript[f, 0]))) (e^(i 2π f t)) dt((|)^2) 
                                              =    (Underscript[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, from the above:
                            PSD(f)=  (Underscript[Lim, T→∞]) ((|A(|^2))/T)((|(∫_0^T) (e^(i 2π (f -(Subscript[f, 0]) ) t)) dt|)^2)
and with                   (Underscript[Lim, T→∞]) (1/T)((|((sin(π T (f-(Subscript[f, 0]))))/(π  (f-(Subscript[f, 0]))))|)^2)  =   δ(f-(Subscript[f, 0]))  ( review delta function if necessary)
Then we obtain:          PSD(f)= (Underscript[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 as expected. 
                        (∫_(-∞)^∞)PSD(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)

6.3 Numerical method to calculate PSD

From the discussion of discrete Fourier transform
                       (∫_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)

Then, based on the definition of  PSD: 
                 PSD[(Subscript[f, p])]= (lim_(T→∞)) (1/T)|(∫_0^T)s(t) (e^(i 2π (Subscript[f, p]) t)) dt  ((|)^2)

The numerical estimate is:   PSD[(Subscript[f, p])]≈ (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.

The unit of the PSD is signal^2/frequency. That’s why it’s call power density. Some example: if s is voltage, |s((|)^2) ≡|V(|^2) the unit of PSD is (V^2)/Hz.  If s is electric current as function of time, it is (A^2)/Hz. Remember that it doesn’t have to be power in the physical sense power=energy/unit time. It can be any quantity that we want to determine its squared amplitude fluctuation vs. time. If it is unit less, such as relative measurement, we just have unitless/Hz.

7. Modulation of harmonic carrier (not for 3340)

8. Advance signal processing techniques (not for 3340)

9. System with nonlinearity (not for 3340)