Sampled time and/or frequency

Discrete time Fourier Transform (DTFT)

The signal $x(n)$ is discrete whereas the resulting frequency spectrum is considered as continuous (arbitrary signals).

• Analysis or direct transform:
    $$X(\omega)=\sum_{n=-\infty}^{\infty} x(n) e^{-j\omega n}$$ (37)

• Synthesis or inverse transform:
  

  $$x(n) = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(\omega) e^{j\omega n} d\omega$$ (38)

  $$= \frac{1}{2\pi} \int_{-0.5}^{0.5} X(f) e^{j 2\pi f n} df$$ (39)

  
  
  note the range here. $f$ is the normalised frequency.

The effect of time domain sampling on the spectrum (Sampling theorem)

Figure 5: Effect of sampling on the spectrum.

What effect has time domain sampling on the frequency spectrum²? Imagine we have an analog spectrum $X_a(F)$ from a continuous signal $x(t)$. We want to know how the spectrum $X(F)$ of the discrete signal $x(n)$ looks like.

  $$X(F) \Leftrightarrow X_a(F)$$ (40)

The signal is represented by $x(n)$ so that we can equate the Fourier transforms of the sampled spectrum $X(F)$ and of the analogue spectrum $X_a(F)$.

  $$\int_{-0.5}^{0.5} \underbrace{X(f)}_{\mbox{sampled}} e^{j 2 \pi f... ...nt_{-\infty}^{+\infty} \underbrace{X_a(F)}_{\mbox{cont}} e^{j 2\pi n F/F_s} dF$$ (41)

Obviously $X(f)$ must be different to accommodate the different integration ranges. The trick is now to divide the integral on the right hand side of Eq. 41 into chunks to make it compatible to the range on the left hand hand side.

Remember that the normalised frequency is $f=F/F_s$ which allows us to change the integration to analogue frequency on both sides:

  $$\frac{1}{F_s} \int_{-F_s/2}^{F_s/2} X(\frac{F}{F_s}) e^{j 2 \pi n F/F_s} dF = \int_{-\infty}^{+\infty} X_a(F) e^{j 2\pi n F/F_s} dF$$ (42)

and now we divide the right hand side into chunks of $F_s$ which corresponds to the integration range on the left hand side.

$$\int_{-\infty}^{+\infty} X_a(F) e^{j 2\pi n \frac{F}{F_s}} dF = \sum_{k=-\infty}^{\infty} \int_{-\frac{1}{2} F_s + k F_s}^{+\frac{1}{2} F_s + k F_s} X_a(F) e^{j2\pi n \frac{F}{F_s}} dF$$ (43)

$$= \sum_{k=-\infty}^{\infty} \int_{-\frac{1}{2} F_s}^{+\frac{1}{2} F_s} X_a(F-k F_s) \underbrace{e^{j2\pi n\frac{F}{F_s}}}_{k F_s \mbox{ omit.}} dF$$ (44)

$$= \int_{-\frac{1}{2} F_s}^{+\frac{1}{2... ...-k F_s)}_{=X(F) \mbox{ of Eq.~\ref{sampl_ana_dig}}} e^{j2\pi n\frac{F}{F_s}} dF$$ (45)

This gives us now an equation for the sampled spectrum:

$$X(F/F_s) = F_s \sum_{k=-\infty}^{\infty} X_a(F-kF_s)$$ (46)

$$X(f) = F_s \sum_{k=-\infty}^{\infty} X_a[(f-k)F_s]$$ (47)

This equation can now be interpreted. In order to get the sampled spectrum $X(F)$ we need to make copies of the analog spectrum $X_a(F)$ and place these copies at multiples of the sampling rate $F_s$ (see Fig. 5). This illustrates also the sampling theorem: if the bandwidth of the spectrum is wider than $F_s/2$ then the copies of the analogue spectrum will overlap and reconstruction would be impossible. This is called aliasing. Note that it is not necessary bad that the spectrum of the analogue signal lies within the range of the so called “base band” $-F/2 \ldots F/2$. It can also lie in another frequency range further up, for example $-F/2 + 34 \ldots F/2 +34$ as long as the bandwidth does not exceed $F_s/2$. If it is placed further up it will automatically show up in the baseband $-F/2 \ldots F/2$ which is called “fold down”. This can be used for our purposes if we want to down mix a signal.

Figure 6: Mapping of frequencies in the analogue domain ($F$) and sampled domain ($f$). $F_s$ is the sampling rate.

With the insight from these equations we can create a plot of how analogue frequencies map onto sampled frequencies. Fig 6 shows how the analogue frequencies $F_\textrm{analogue}$ map on the normalised frequencies $f_\textrm{sampled}$. As long as the analogue frequencies are below $F_s/2$ the mapping is as usual as shown in Fig 6A. Between $F_s/2$ and $F_s$ we have an inverse mapping: an increase in analogue frequency causes a decrease in frequencies. Then, from $F_s$ we have again an increase in frequencies starting from DC. So, in general if we keep a bandlimited signal within one of these slopes (for example from $F_s \ldots F_s + 1/2 F_s$ as shown in Fig 6B) then we can reproduce the signal.

This leads us to the generalised Nyquist theorem: if a bandpass filtered signal has a bandwidth of $B$ then the minimum sampling frequency is $F_s = 2B$.

Discrete Fourier Transform (DFT)

So far we have only sampled in the time domain. However, on a digital computer the Fourier spectrum will always be a discrete spectrum.

The discrete Fourier Transform (DFT) is defined as:

  $$X(k) = \sum_{n=0}^{N-1} x(n) e^{-j2\pi kn / N} \qquad k=0,1,2, \ldots, N-1$$ (48)

where $N$ is the number of samples in both the time and frequency domain.

The inverse discrete Fourier Transform (IDFT) is defined as:

  $$x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X(k) e^{j 2\pi kn/N} \qquad n=0,1,2,\ldots, N-1$$ (49)

Figure 7: Effect of sampling in the frequency domain. The inverse in the time domain $x(n)$ contains repetitions of the original signal.

What is the effect in the timedomain of this discretisation³? We start with the continuous Fourier transform and discretise it into N samples in the frequency domain:

  $$X\left(\frac{2\pi}{N}k\right) = \sum_{n = -\infty}^{\infty} x(n) e^{-j\frac{2\pi}{N}kn} \qquad k = 0, \ldots, N-1$$ (50)

Let's subdivide the sum into chunks of length $N$:

$$X\left(\frac{2\pi}{N}k\right) = \sum_{l = -\infty}^{\infty} \sum_{n=lN}^{lN+N-1} x(n)e^{-j \frac{2\pi}{N}kn}$$ (51)

$$= \sum_{l=-\infty}^{\infty} \sum_{n=0}^{N-1} x(n-lN) e^{-j \frac{2\pi}{N} kn}$$ (52)

$$= \sum_{n=0}^{N-1} \underbrace{\sum_{l=-\infty}^{\infty} x(n - lN)}_{\mbox{Periodic repetition!}} e^{-j \frac{2\pi}{N}kn}$$ (53)

We note the following:

• Ambiguity in the time domain
• The signal is repeated every N samples

Practically this repetition won't show up because the number of samples is limited to $N$ in the inverse transform. However, for operations which shift signals in the frequency domain it is important to remember that we shift a periodic time series. If we shift it out at the end of the array we will get it back at the start of the array.

Figure 8: Example of a DFT. The sampled signal $x(n)$ is converted to the spectrum $X(k)$. Both have $N$ samples. Note the redundancy: the spectrum is mirrored around $N/2$. The first value $X(0)$ is the DC value of the signal. However this is the only sample which is not mirrored.

Fig. 8 shows an example of a DFT. It is important to know that the spectrum is mirrored around $N/2$. DC is represented by $X(0)$ and the Nyquist frequency $F_s/2$ is represented by $X(N/2)$. The mirroring occurs because the input signal $x(n)$ is real. This is important if one wants to modify the spectrum $X(F)$ by hand, for example to eliminate 50Hz noise. One needs to zero two elements of $X(F)$ to zero. This is illustrated in this python code:

import scipy as sp
yf=sp.fft(y)
# the sampling rate is 1kHz. We've got 2000 samples.
# midpoint at the ifft is 1000 which corresponds to 500Hz
# So, 100 corresponds to 50Hz
yf[99:101+1]=0;
# and the mirror
yf[1899:1901+1]=0;
#
yi=sp.ifft(yf);

This filters out the 50Hz hum from the signal $y$ with sampling rate 1000Hz. The signal yi should be real valued again or contain only very small complex numbers due to numerical errors.

Properties of the DFT

• Periodicity:
  

  $$x(n + N) = x(n)$$ (54)

  $$X(k + N) = X(k)$$ (55)

  
  

• Symmetry: if $x(n)$ real:
    $$x(n)\mbox{ is real} \Leftrightarrow X^*(k) = X(-k) = X(N-k)$$ (56)
  This is important when manipulating $X(k)$ by hand.

• Time Reversal:
  

  $$x(n) \leftrightarrow X(k)$$ (57)

  $$x(-n) \leftrightarrow X(N - k)$$ (58)

  
  

• Circular convolution:
    $$X_{1}(k) X_{2}(k) \leftrightarrow x_1(n) * x_2(n)$$ (59)
  with
    $$x_{3}(m) = \sum_{n = 0}^{N - 1} x_{1}(n) x_{2}(m - n)$$ (60)

More useful equations are at Proakis and Manolakis (1996, pp.415).

Problems with finite length DFTs

Figure 9: The effect of windowing on the DFT.

  $$x_w(n) = x(n) \cdot w(n)$$ (61)

where $x$ is the input signal and $w(n)$ represents a function which is $1$ from $n=0$ to $n=L-1$ so that we have $L$ samples from the signal $x(n)$.

To illustrate the effect of this finite sequence we introduce a sine wave $x(n) = \cos \omega_0 n$ which has just two peaks in a proper Fourier spectrum at $-\omega$ and $+\omega$.

The spectrum of the rectangular window with the width $L$ is:

  $$W(\omega) = \frac{sin(\omega L /2)}{sin (\omega / 2)} e^{-j\omega (L -1)/2}$$ (62)

The resulting spectrum is then:

  $$X_{3}(\omega) = X(\omega) * W(\omega)$$ (63)

Because the spectrum of $X(\omega)$ consists of just two delta functions the spectrum $X_{3}(\omega)$ contains the window spectrum $W(\omega)$ twice at $-\omega$ and $+\omega$ (see Fig. 9). This is called leakage. Solutions to solve the leakage problem? Use Windows with a narrower spectrum and with less ripples (see FIR filters).

Fast Fourier Transform

We can rewrite the DFT (Eq. 48) in a slightly more compact form:

  $$X(k) = \sum_{n = 0}^{N-1} x(n) W_{N}^{kn}$$ (64)

with the constant:

  $$W_{N} = e^{-j 2\pi/N}$$ (65)

The problem with the DFT is that it needs $N^2$ multiplications. How can we reduce the number of multiplications? Idea: Let's divide the DFT in an odd and an even sequence:

$$x(2m)$$ (66)

$$x(2m + 1), \qquad m = 0, ..... , \frac{N}{2} - 1$$ (67)

which gives us with the trick $W_N^{2mk} = W_{N/2}^{mk}$ because of the definition Eq. 65.

$$X(k) = \sum_{m=0}^{N/2-1} x(2m) W_{N}^{2mk} + \sum_{m=0}^{N/2 - 1} x(2m + 1) W_{N}^{k (2m + 1)}$$ (68)

$$= \sum_{m=0}^{N/2-1} x(2m) W_{N/2}^{mk} + W_{N}^{k} \sum_{m = 0}^{N/2 - 1} x(2m + 1) W_{N/2}^{mk}$$ (69)

$$= F_{e}(k) + W_{N}^{k} F_{o}(k)$$ (70)

$F_e$ and $F_o$ have both half the length of the original sequence and need only $(N/2)^2$ multiplication, so in total $2 \cdot (N/2)^2 = \frac{N^2}{2}$. Basically by dividing the sequence in even and odd parts we can reduce the number of multiplications by 2. Obviously, the next step is then to subdivide the two sequences $F_{e}(k)$ and $F_{o}(k)$ even further into something like $F_{ee}(k), F_{eo}(k), F_{oe}(k)$ and $F_{oo}(k)$.

Figure 10: Illustration of the Fast Fourier Algorithm. The sequence of $N$ samples is recursively divided into subsequences of odd and even samples.

In general the recipe for the calculation of the FFT is:

  $$X_{i}(k) = X_{ie}(k) + W_{L}^{k} X_{io}(k)$$ (71)

$W_{L}^{k}$ is the phase factor in front of the odd sequence. This is continued until we have only two point ($N=2$) DFTs (see Eq. 64):

$$\mbox{DC:} \qquad X(0) = x(0) + \underbrace{W_2^0}_{1} x(1) = x(0) + x(1)$$ (72)

$$\mbox{Nyquist frequ.:} \qquad X(1) = x(0) + \underbrace{W_2^1}_{-1} x(1) = x(0) - x(1)$$ (73)

A two point DFT operates only on two samples which can represent only two frequencies: DC and the Nyquist frequency which makes the calculation trivial. Eq. 72 is an averager and Eq. 73 is basically a differentiator which gives the max output for the sequence $1,-1,1,-1,\ldots$. Fig. 10 illustrates how to divide the initial sequence to arrive at 2 point DFTs. In other words the calculation of the full DFT is done by first calculating $N/2$ 2 point DFTs and recombining the results with the help of Eq. 71. This is sometimes called the “Butterfly” algorithm because the data flow diagram can be drawn as a butterfly. The number of complex multiplications reduces in this approach to $N \log_2 N$ which is actually the worst case scenario because many $W_{L}^{k}$ usually turn out to be $1,-1,j,-j$ which are just sign inversions or swaps of real and imaginary parts. A clever implementation of this algorithm will be even faster.

In summary the idea behind the FFT algorithms is to divide the sequence into subsequences. Here we have presented the most popular radix 2 algorithm. The radix 4 is even more efficient and there are also algorithms for divisions into prime numbers and other rather exotic divisions. However, the main idea is always the same: subsample the data in a clever way so that the final DFT becomes trivial.

github / contact