Numerical Methods for Fourier Analysis of LTI
- Part 1

Discrete Fourier Transform

ECE 3336/3337/3340/4375

Han Q. Le (c)        Also for ECE6323

1 Introduction review of LTI (see ppt file)

Time domain

Frequency domain

The harmonic transfer function is related to the Laplace transfer function:
                              H[ω]= G[i ω]=|G[i ω]| (e^(i Φ[G[i ω]]))                            (1.1.1)

- Harmonic transfer function: assume no beginning: in practice, signal at long time steady state.
- Very useful to study system nonlinearity by looking at higher-order harmonics or combinations of multiple input harmonics
- Key terms here: amplitude (magnitude) of the response: |G[i ω]| and phase term of the response: (e^(i Φ[G[i ω]])); Φ[G[i ω]] is called the phase.
- If amplitude is >1: amplification effect. If If amplitude is <1: attenuation effect. Its frequency-dependence behavior is essential; e. g. filter.
- Phase can be advanced or retarded (delay). Let input:
            (Subscript[x, i])[t]=Cos[ω t]     and output:    (Subscript[x, o])[t]=Cos[ω t+ϕ] =Cos[ω (t+(ϕ/ω))]    
we see that the term (ϕ/ω) acts like the time for signal translation, which can be advanced or delay. This is corresponding to the lag-lead concept in transient response that we learned in the previous chapter.

Manipulate[Plot[{Cos[2 π t], Cos[2 π (t+τ)]}, {t,0,5},PlotStyle-> {{Thick, Red},{Thick,Blue}}],{{τ,0},- 1, 1}]

Notice the concept of which signal that comes first (advanced) relative with which signal that is delayed.

2. Basic DFT

2.1 Review and general introduction

We have seen this in previous lectures on Fourier transform:

If we know, or have a good model of the signal of interest, analytic method can give us a formula of the Fourier transform. But as we discussed in the lecture on power spectral density, the real interest of Fourier analysis is for empirical application, where the signal property is not completely known. In fact, numerical Fourier method is applied to data in order to determine the characteristic properties of a system so that it can be modeled. In other words, empirical analysis usually comes first, and modeling comes after.

So, how is the data of a signal or system obtained? By discrete sampling. For example, below are data points of a signal.

. What is the true signal like? Here it is:

As we can see, it is not always obvious to us what the signal really is based on measurements. Mathematical methods must be applied to obtained the best inference of the true signal or system properties from the finite and limited data.

2.2 Concept of numerical Fourier transform

For signals and systems with periodic properties or harmonic responses, we know that Fourier transform is very useful. But we only have data points, how do we obtain the best estimate of its Fourier transform?

More precisely, we can rephrase the problem like this: If we have a finite series of measurement (or sampled points): (Subscript[{(Subscript[s, k])}, k=1]^M) for signal s[t], how do we find its Fourier transform?
                                      (∫_(-∞)^∞)s(t)(e^(i ω t))dt   (or  (∫_(-∞)^∞)s(t)(e^(i 2 π f t))dt  ) ?                 (2.2.1)
(Note: the usual factor (1/(2 π)) or (1/(Sqrt[2 π])) are just by convention and can be put in later. They are removed for now to allow clearer presentation of the approach).

Naturally, it makes sense that we can substitute the integration with a finite summation as an estimate.  Let say the series (Subscript[{(Subscript[s, k])}, k=1]^M) are measured with regular interval (either time, space, or whatever), then the series (Subscript[{(Subscript[s, k])}, k=1]^M) actually represent pairs of numbers:  
                                                    (Subscript[{((Subscript[t, k]) , (Subscript[s, k]))  }, k=1]^M)  ≡(Subscript[{( (k-1)Δt, (Subscript[s, k]))  }, k=1]^(M  ))              (2.2.2)
where Δt is the sampling interval. We can take a numerical estimate of the Fourier transform of this data  by approximating 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))            (2.2.3)

Fortunately, if f is divisible with respect to (1/(M Δt)), i. e., it is an integral multiple of that quantity, then there is a very efficient algorithm to perform the summation of (2.2.3), in fact it is called FFT or Fast Fourier Transform for that reason. Then, it appears this problem is solved. However, it is not so simple. Consider this example:

If we apply (2.2.3), how certain are we about the estimate of its FT? Below is the actual signal:

The estimate, if we apply (2.2.3) will be totally off. The fact is, variable f in  (2.2.3) can not really be treated as a continuous variable over an infinite range. Hence, what important in learning and using DFT is to know the limit of its validity and how to apply it appropriately.

For the above example, it is clear that if the sample rate is too slow compared to the signal oscillation, we cannot be certain about the estimated FT for f beyond certain range. In more specific terms, by a theorem known as Nyquist’s theorem, we can’t tell if there is any frequency higher than half of the sampling frequency. In other words, if we sample a signal at a rate of say 1,000,000 samples per second (or 1 MHz), we can’t tell if there is a frequency at 0.5 MHz or higher. It takes at least 2 samples per cycle to detect the frequency.

Likewise, no sampling episode can be infinite. The longest time we measure is  M Δt. This means if there is any trend that varies slowlier than the episode M Δt, we can’t tell it either.

2.3 Discrete Fourier transform

Thus, the frequencies that can be valid for a FT estimate have only have a finite range from smallest to largest. Futhermore, to obtain its numerical values, we also need an array of discretized  f  where the FT is computed. Let this array be a series of sample frequencies: (p-1) Δf where p is an integer and Δf is the frequency incremental interval. It turns out that the most efficient and fastest algorithm to compute (2.2.3) is to choose:
                                                                Δf =(1/(MΔt))                   (2.3.1)
 Then, (2.2.3) becomes:
                                          (∫_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)) (2.3.2)
and this is known as discrete Fourier transform (DFT).

In the fast Fourier transform (FFT) numerical technique, quantity (Subscript[s, k])(e^(i 2π (p-1) (k-1) /M)) can be computed with a highly efficient algorithm. Generally, DFT (Subscript[v, s]) of a list (Subscript[u, r]) of length n is defined to be:
                                                             (Subscript[v, s])=(1/(Sqrt[n]))(∑_(r=1)^n)(Subscript[u, r])(e^(2π i (r-1)(s-1)/n))
For Mathematica, this is simply done with function Fourier[].
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)

Do not forget that Mathematica uses (1/(Sqrt[2π])) default convention for its FourierTransform[], hence, if both are used somehow in a calculation (when both data and analytic model are involved), for consistency, the above formula should have (1/(Sqrt[2π])) on both sides (and  ω=2 π (Subscript[f, p]))
                                 (1/(Sqrt[2π]))(∫_0^(M Δt))s(t)(e^(i 2π  (Subscript[f, p]) t))dt  ≈ ((Δt)/(Sqrt[2π]))  (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π (p-1) (k-1) /M))
                                                         = Δt  (Sqrt[(M/(2 π))]) FFT((Subscript[{(Subscript[s, k])}, k=1]^M);p)

Most commercial computational software can perform this operation very fast and efficient. We will see examples with Mathematica in subsequent section. But first, let’s look at few basic properties to familiarize ourselves with numerical DFT.

2.4 Some basic properties

2.4.1 A simple example

Consider a pure harmonic signal Sin[2 π f t], let’s sample it and find its numerical Fourier transform:

Demo: basic illustration of harmonic signal

Notice something? In the DFT output, the peaks correspond to the frequency of the signal for sure, but why there are 2 peaks that are symmetric across the middle point of the array?
The reason is this. Consider the Fourier component S[(Subscript[f, p])] of frequency (Subscript[f, p])=(p-1) Δf :
                                                    S[(Subscript[f, p])]=Δt  (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π (p-1) (k-1) /M))
Let’s substitute M-p+2  for p:
                                             S[(Subscript[f, M-p+2])]=Δt  (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π (M-p+2-1) (k-1) /M))    
                                                              =Δt  (∑_(k=1)^M)(Subscript[s, k]) (e^(i 2π (M-p+1) (k-1) /M)) 
                                                              =Δt  (∑_(k=1)^M)(Subscript[s, k])(e^(i 2π M (k-1) /M))(e^(-i 2π (p-1) (k-1) /M))
                                                              =Δt  (∑_(k=1)^M)(Subscript[s, k])(e^(-i 2π (p-1) (k-1) /M)) ~(∫_0^(M Δt))s(t)(e^(-i 2π  (Subscript[f, p]) t))dt

(∑_(k=1)^M)(Subscript[s, k])(e^(i 2π (p-1) (k-1) /M))

In other words, S[(Subscript[f, M-p+2])] is actually the negative frequency component corresponding to (Subscript[f, p]):
                                                    S[(Subscript[f, M-p+2])]= S[-(Subscript[f, p])]   
Note that p=1 is just the zero frequency component, and (Subscript[f, M-p+2])=(Subscript[f, M+1])  does not exist. This should be expected, since there is only one zero-frequency component. 
Then                                              S[(Subscript[f, M])]= S[-(Subscript[f, 2])]

2.4.2 A simple example - continue

Hence, if we cut the second half of the output DFT and put on the first half, and associate the array with frequencies on both minus and plus side:            
                       ({{f, -> , -(M/2)Δf, -((M/2)-1)Δf, ..., -2 Δf,  - Δf, 0,  Δf, ..., ((M/2)-2)Δf, ((M/2)-1)Δf}, {p, -> , (M/2)+1, (M/2)+2, ..., M-1, M, 1, 2, .., (M/2)-1, (M/2)}})
We notice that if M is even, there is a slight asymmetry: one component has to be missed on the extreme positive or extreme negative side. However, this is not a problem really. If it is odd, we can have even component both sides.
If it is even:                          S[(Subscript[f, M-((M/2)+1)+2])]= S[-(Subscript[f, (M/2)+1])]  
                                                S[(Subscript[f, (M/2)+1])]= S[-(Subscript[f, (M/2)+1])]  
What this means is that we can add an extra component on the right to complement with the left most component.

Demo: basic illustration of how to rotate a DFT array to obtain negative frequency part

Demo: basic illustration of harmonic signal spectrum

2.5 Positive frequency, FT and Inverse FT

For real physical time signal, we don’t have negative frequency. But for spatial variables, we can use all spatial frequency in 3D, what if we want to evaluate for negative value of f?
                                                                   (∫_0^(M Δt))s(t)(e^(-i 2π f t))dt   
This is known as inverse Fourier, and 
and inverse DFT gives us:
                                                              (Subscript[w, s])=(1/(Sqrt[n]))(∑_(r=1)^n)(Subscript[u, r])(e^(-2π i (r-1)(s-1)/n))

2.5.1 Inverse DFT demo

Demo: basic illustration of harmonic signal

2.6 Examples of common pulse shape

Consider the signal in the app below:

Demo: FFT of common pulse shape

3. (not for ECE 3340)

4. LTI and signal processing - part 1

4.1

Now we can apply to signal processing. Consider, for example a very simple RC filter:
                                             H[f]=(1/(1+i  2 π f  τ))        or     H[f]=((i  2 π f  τ)/(1+i  2 π f  τ)) 
where                                                   τ=R C

We can also have a typical RLC resonance filter:

H[f]=((((f/(Subscript[f, 0])))^2)/((((f/(Subscript[f, 0])))^2)-1-i  2 π f  τ)) 

 Let’s put a pulse through the filter, we can see what we get.

Demo of RCL filters

5. Signal with carrier (optional for ECE 3340)

5.1 Discussion

Consider a signal s[t] with a carrier (e^(-i (Subscript[ω, 0]) t))s[t], recall from FT theorem, its transform is: 
                                     (∫_(-∞)^∞)(e^(i ω t))(e^(-i (Subscript[ω, 0]) t))s[t]dt =(∫_(-∞)^∞)(e^(i  (ω-(Subscript[ω, 0])) t))s[t]dt =𝔉[s;(ω-(Subscript[ω, 0]))]≡ S[ω-(Subscript[ω, 0])] 
If we don’t know the signal carrier, as illustrated in demo with sound signal, then we have to sample with a sufficiently fine interval to determine the carrier frequency. However, for simulation for example, if we do know the carrier, then the theorem tells us that we have only to do the transform for the signal envelop. Then we can simply shift that transform to the origin at (Subscript[ω, 0])

For real signal, such as Cos[(Subscript[ω, 0]) t]
   (1/2)(∫_(-∞)^∞)(e^(i ω t))((e^(i (Subscript[ω, 0]) t))+(e^(-i (Subscript[ω, 0]) t)))s[t]dt =(1/2) (S[ω-(Subscript[ω, 0])]+ S[ω+(Subscript[ω, 0])] )
   Or:                       (1/2) (S[ (f-(Subscript[f, 0]))]+ S[(f+(Subscript[f, 0]))] )

Numerically:      (Subscript[S, C])[m Δf+(Subscript[f, 0])]=(1/2) (S[ m Δf]+ S[(m Δf+2(Subscript[f, 0]))] )
                           (Subscript[S, C])[m Δf-(Subscript[f, 0])]=(1/2) (S[ m Δf-2 (Subscript[f, 0])]+ S[m Δf] )

However, we may not have S[(m Δf+2(Subscript[f, 0]))] and S[ m Δf-2 (Subscript[f, 0])] because they are outside sampling range. This is illustrated as below

However, a system response function may have complicated dependence on
                                      (Subscript[S, out])[ω]=H[ω]S[ω-(Subscript[ω, 0])]  
                                      (Subscript[S, out])[ω]=H[(Subscript[ω, 0])+Ω]S[Ω]  
                   (Subscript[s, out])[t]=(1/(2 π))(∫_(-∞)^∞)(e^(-i ω t))(Subscript[S, out])[ω] dω=(1/(2 π))(∫_(-∞)^∞)(e^(-i ω t))H[(Subscript[ω, 0])+Ω]S[Ω]  dω   
                             =(1/(2 π))(e^(-i (Subscript[ω, 0]) t))(∫_(-∞)^∞)(e^(-i (ω -(Subscript[ω, 0]))t))H[(Subscript[ω, 0])+Ω]S[Ω]  dω   
                             =(1/(2 π))(e^(-i (Subscript[ω, 0]) t))(∫_(-∞)^∞)(e^(-i Ω t))H[(Subscript[ω, 0])+Ω]S[Ω]  dΩ

Very often, we can expand H[(Subscript[ω, 0])+Ω] as follow:
                  H[(Subscript[ω, 0])+Ω]=|H[(Subscript[ω, 0])+Ω]|(e^(i Φ[(Subscript[ω, 0])+Ω]))
                                    =|H[(Subscript[ω, 0])+Ω]|(e^(i Φ[(Subscript[ω, 0])])) (e^(i Ω (Subscript[τ, 0])))(e^(i (d/2)(Ω^2)))(e^(i χ[Ω]))
 where:                        (Subscript[(Subscript[τ, 0])=((dΦ)/(dω))|, (Subscript[ω, 0])])  ,  (Subscript[d=(((d^2)Φ)/(d(ω^2)))|, (Subscript[ω, 0])])

We will see that each term has a meaning.

5.2 Envelope or group delay

(Subscript[s, out])[t]=(1/(2 π))(e^(-i (Subscript[ω, 0]) (t-((Φ[(Subscript[ω, 0])])/(Subscript[ω, 0])))))(∫_(-∞)^∞)(e^(-i Ω (t-(Subscript[τ, 0]))))|H[(Subscript[ω, 0])+Ω]|(e^(i (d/2)(Ω^2)))(e^(i χ[Ω]))S[Ω]  dΩ  
If d is negligible and to the first order approximation, |H[(Subscript[ω, 0])+Ω]| is nearly constant, then:
               (Subscript[s, out])[t]=(1/(2 π))(e^(-i (Subscript[ω, 0]) (t-((Φ[(Subscript[ω, 0])])/(Subscript[ω, 0])))))|H[(Subscript[ω, 0])]|(∫_(-∞)^∞)(e^(-i Ω (t-(Subscript[τ, 0]))))S[Ω]  dΩ   
              (Subscript[s, out])[t]=(1/(2 π))(e^(-i (Subscript[ω, 0]) (t-((Φ[(Subscript[ω, 0])])/(Subscript[ω, 0])))))|H[(Subscript[ω, 0])]|(Subscript[s, in])[t-(Subscript[τ, 0])]      
The meaning is clear, ((Φ[(Subscript[ω, 0])])/(Subscript[ω, 0])) is the phase delay term, and (Subscript[τ, 0]) is the envelop delay term (if for traveling wave, it is also known as group delay term). It is rarely ((Φ[(Subscript[ω, 0])])/(Subscript[ω, 0]))=(Subscript[τ, 0]) , which means that the system introduces a relative shift between the carrier and the envelope as illustrated below.
This is the case when we have no need for any computation because the output is obtained directly from the input.

5.3 Quadratic dispersion term - pulse distortion

For the next term, it may result in

(∫_(-∞)^∞)(E^(I Ω t)) (E^((-(t^2))/(2(T^2))))\[DifferentialD]t

(∫_(-∞)^∞)(E^(-I Ω τ)) (E^(((-1)/2)(Ω^2) (P^2)))\[DifferentialD]Ω

(((E^(-((τ^2)/(2 (P^2))))) (Sqrt[2 π]))/(Sqrt[(P^2)])) /. (P^2)-> (T^2)+I (d/2) /. τ-> t-(Subscript[τ, 0])

Hence, the output is:      
            (Subscript[s, out])[t]=|H[(Subscript[ω, 0])]|(e^(-i (Subscript[ω, 0]) (t-((Φ[(Subscript[ω, 0])])/(Subscript[ω, 0])))))(Sqrt[(1/(1+((i d)/(2 (T^2)))))])(e^(-((((t-(Subscript[τ, 0])))^2)/(2 ((T^2)+i (d/2))))))

5.4 General computation

When other higher order terms are not neglible, we have to do numerical calculation:
                    (Subscript[s, out])[t]=(1/(2 π))(e^(-i (Subscript[ω, 0]) (t-((Φ[(Subscript[ω, 0])])/(Subscript[ω, 0])))))(∫_(-∞)^∞)(e^(-i Ω (t-(Subscript[τ, 0]))))|H[(Subscript[ω, 0])+Ω]|(e^(i (d/2)(Ω^2)))(e^(i χ[Ω]))S[Ω]  dΩ       
                             =(e^(-i (Subscript[ω, 0]) (t-((Φ[(Subscript[ω, 0])])/(Subscript[ω, 0])))))(𝔉^(-1))[U;(t-(Subscript[τ, 0]))](∫_(-∞)^∞)(e^(-i Ω (t-(Subscript[τ, 0]))))|H[(Subscript[ω, 0])+Ω]|(e^(i (d/2)(Ω^2)))(e^(i χ[Ω]))S[Ω]  dΩ     
  where we define                U[Ω]≡|H[(Subscript[ω, 0])+Ω]|(e^(i (d/2)(Ω^2)))(e^(i χ[Ω]))S[Ω] ≡ G[Ω] S[Ω]   
                                                         G[Ω]≡|H[(Subscript[ω, 0])+Ω]|(e^(i (d/2)(Ω^2)))(e^(i χ[Ω]))

We note that usually (Subscript[ω, 0])>>(Subscript[Ω, M]), which is ~ the bandwidth of envelop S[Ω], and hence unless |H[(Subscript[ω, 0])+Ω]| is a highly varying function vs (Subscript[ω, 0]), which is extremely unlikely to be the case, G[Ω] is really a low-frequency function (i. e. envelop bandwidth) and it means that for computation purpose, there is no need for sampling at the carrier frequency but only sufficient for the envelop.