Lab VII The Discrete Fourier Transform

As emphasized in the previous labs, sinusoids are an important part of signal analysis, since many signals that occur in the real world are composed of sinusoids. For example, many musical signals can be approximately described as sums of sinusoids, as can some speech sounds (vowels in particular). It turns out that any periodic signal can be written exactly as a sum of amplitude-scaled and phase-shifted sinusoids. Equivalently, we can use Euler’s inverse formulas to write periodic signals as sums of complex exponentials, a mathematically more convenient description, and the one that we will adopt in this course. The description of a signal as a sum of sinusoids or complex exponentials is known as the spectrum of the signal.

So far, we have typically thought of signals as time-varying quantities, like x(t) or x[n]. When we plot these signals, we generally place time along the horizontal axis and signal value along the vertical axis. The idea behind the frequency-domain representation of a signal is similar. Rather than plotting signal value versus time, we plot a spectral value versus frequency. Doing this involves a transformation of the signal. The figure below shows an example of a time-domain and frequency-domain representation of a signal (in fact the modulus and phase of the coefficients $X_k$ of the Fourier series decomposition of the signal). Note that we can think of the result of the transform as a signal as well, a signal whose independent variable is frequency rather than time.

Why do we need another representation for a signal ? Isn’t the usual time-domain representation enough ? It turns out that spectral (or frequency-domain) representations of signals have many important properties. First, a frequency domain representation may be simpler than a time-domain representation, especially in cases where we cannot write an analytic expression for the signal. Second, a frequency-domain representation of a signal can often tell us things about the signal that we would not know from just the time-domain signal. Third, a signal’s spectrum provides a simple way to describe the effect of certain systems (like filters) on that signal. There are many more uses for frequency-domain representations of a signal, and we will examine some of them throughout this course. Spectral representations are one of the most central ideas in signals and systems theory, and can also be one of the trickiest to understand.

Consider the following problem. Suppose that we have a signal that is actually the sum of two different signals. Further, suppose that we would like to separate one signal from the other, but the signals overlap in time. If the signals have frequency-domain representations that do not overlap, it is still possible to separate the two signals. In this way, you can see that frequency-domain representations provide another “dimension” to our understanding of signals.

In brief, the main focus of this lab is to learn how to compute the spectrum of (discrete-time) signals.

VII.1 Using the DFT to evaluate the coefficients of the Fourier series development of a periodic or time-limited signal

Background

In order to evaluate the coefficients of the Fourier series development of a continuous-time (finite) power signal x(t), supposed to be periodic of fundamental period T, or time-limited with a finite duration T, we have to follow these steps :

1. Sampling x(t) with a sampling period $t_e=\Delta t$ : $$x\left( {k\Delta t} \right) = x\left( t \right)\delta _{\Delta t} \left( t \right) \equiv x\left[ k \right],\;k \in \mathbb{Z}$$

2. Windowing x[k] to delimit a finite number of samples of the initial (infinite) sequence, to keep N samples which define an overall sampling duration $T_0=N \Delta t$. The ideal windowing takes place on an integer number of the fundamental period, so that : $$T_0 = N\Delta t = mT$$

3. Discretization of the frequency axis : We start from the expression of the SF coefficient $X_k$. Although $X_k$ is usually defined and calculated on a single period T of x(t), there is no problem to use an integer multiple of T : $$X_k = \frac{1}{T}\int\limits_{t_1 }^{t_1 + T} {x(t)e^{ - j2\pi \frac{k}{T}t} dt} = \frac{1}{{mT}}\int\limits_{t_1 }^{t_1 + mT} {x(t)e^{ - j2\pi \frac{k}{{mT}}t} dt}$$ While it seems that the frequency of the fundamental component will be reduced to $f'_0=\frac{f_0}{m}$, it is easy to verify that the additional points are all zeros, such that the fundamental component will still appear at the right frequency $f_0$ !

4. Discretize the expression of the $X_k$ coefficients : ParseError: KaTeX parse error: Undefined control sequence: \displaylines at position 2: \̲d̲i̲s̲p̲l̲a̲y̲l̲i̲n̲e̲s̲{ X_W \left[ …

One can thus obtain the exact value of the SF development : $$\displaystyle{X_k = X_W \left[ k \right] = \frac{1}{N}DFT\left( {x\left[ n \right]} \right)}$$ Provided that :

1. The sampling of x(t) respects the Shannon-Nyquist condition for a sufficiently fast sampling, defined for a (low-pass) band-limited signal of bandwidth B Hz as $f_s = \frac{1}{{\Delta t}} \ge 2B$.

2. The truncation has been adequately done on an integer number of periods : $T_0 = N\Delta t = mT$.

This way, it is generally possible to evaluate the exact value of the SF development coefficients using the Fast Fourier Transform (FFT) on any computer. Using Octave/MatLab, this algorithm is launched using the instruction fft, which realizes the DFT defined as : $$X_{k + 1} = \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {x[n + 1]e^{ - j2\pi \frac{k}{N}n} } \quad k,n = 0,...,N - 1$$

Remarks

• The first element of a vector is indexed by 1 in Octave/MatLab, and there is no negative values for vector or matrix indexes.

• Don’t forget to divide by the total number of samples N to find the correct amplitude of the SF coefficients.

• The frequency interval between two neighbouring points, which is $\Delta f = \frac{1}{T_0}$, is not used, nor computed by the fft. You should explicitly generate a frequency support vector if you want to plot the spectrum vs the right frequency.

• The coefficients which correspond to negative frequency indexes k are located in the second half of the output vector. In order to obtain the standard bilateral representation centered on the zero frequency, the raw DFT output vector should be rearranged, using e.g. the dedicated function fftshift.

• The FFT algorithm is particularly efficient when N is a power of 2, a little bit less if N is even. It is rather slow when N is a prime. You should particularly see this difference using 2D DFT.

Exercise VII.1

1. Plot in the first panel of a multi-graph window the representation of the periodic (power) signal $s_1(t) = 1+2 sin(2\pi f_1 t)$, using $f_1 = 3 Hz$ and a total duration $T_0=1s$. Use the DFT to evaluate the complex coefficients of the FS decomposition, using a total sampling duration of exactly one period of the discrete sinusoid (hint : pay attention to the fact that the last point of the vector should not be equal to the first ! ). Plot their real and imaginary part in the second and third graph, and their modulus and phase in a third and fourth graph. The four graphs should be correctly indexed horizontally with the right frequencies. Verify that only three coefficients are different from zero, and if not, try to understand why ( did you really sampled on one or an integer number of periods ? ). Show your results to your instructor and discuss them with him.

2. Redo the previous multi-graph window and evaluate the SF coefficients $X_k$ of the same periodic signal using a total sampling duration equals to an integer multiple m≠1 of the fundamental period. Explain what happens on the spectrum.

3. Redo the previous multi-graph window and evaluate the $X_k$ of the same periodic signal using a total sampling duration equals to (m+0.25) times the fundamental period. Explain what happens on the spectrum.

4. For which value of $T_0$ (in terms of the integer m) do you expect the apodization effect to be the most prominent ?

VII.2 Using the DFT to evaluate the Fourier Transform of a signal

Background

In order to evaluate the Fourier transform of a continuous-time signal x(t), which may have an infinite duration, we have to follow these steps :

1. Sampling x(t) with a sampling period $t_e=\Delta t$ : $$x\left( {k\Delta t} \right) = x\left( t \right)\delta _{\Delta t} \left( t \right) \equiv x\left[ k \right],\;k \in \mathbb{Z}$$

2. Windowing x[k] to delimit a finite number of samples of the initial sequence, to keep N samples which covers the total signal duration $T_0=N \Delta t$.

3. Discretization of the frequency axis : Knowing that the Fourier Transform $X(f)$ of a discrete-time signal is periodic with a fundamental period $f_e = \frac{1}{\Delta t}$, it is enough to sample $X_W(f)$ on a single period using N points (N input values, N output values !).

   Therefore, the frequency interval $\Delta f$ is automatically determined : $$N\Delta f = \frac{1}{{\Delta t}}\quad \to \quad \Delta f = \frac{1}{{N\Delta t}} = \frac{1}{{T_0 }}$$ The Fourier transform $X(f) = \int\limits_{ - \infty }^{ + \infty } {x(t)e^{ - j2\pi ft} dt}$ can now be discretized and approximated as follows : ParseError: KaTeX parse error: Undefined control sequence: \displaylines at position 2: \̲d̲i̲s̲p̲l̲a̲y̲l̲i̲n̲e̲s̲{ \mathop {\l… which gives : $$\mathop {\left. {X(f)} \right|}\nolimits_{f = \frac{n}{{T_0 }}} \approx \mathop {\left. {X_W (f)} \right|}\nolimits_{f = \frac{n}{{T_0 }}} \approx \frac{T_0}{N} DFT\left( {x\left[ k \right]} \right)$$

Remarks

1. The sampling of x(t) respects the Shannon-Nyquist condition for a sufficiently fast sampling, defined for a (low-pass) band-limited signal of bandwidth B Hz as $f_s = \frac{1}{{\Delta t}} \ge 2B$.

2. Frequency interval : The frequency interval $\Delta f$ between two neighboring points, which is given by $\Delta f = \frac{1}{{N\Delta t}} = \frac{1}{{T_0 }}$ should not be confused with the sampling frequency $f_s$, since $N\Delta f = \frac{1}{{\Delta t}} = f_s \quad \to \quad \Delta f = \frac{1}{{N\Delta t}} = \frac{{f_s }}{N}$.

   An important rule : Since $\Delta f = \frac{1}{T_0}$, the longest is $T_0$ , or the highest is N since $T_0 = N \Delta t$, the smallest $\Delta f$ will be, and the better the frequency « resolution » of the DFT result, knowing that a high resolution will allow to distinguish e.g., two adjacent peaks in the spectrum. Sometimes, one may add a sequence of zeros to artificially increase the total duration of the input temporal sequence, a procedure called the “zero-padding” trick.

3. Frequency scaling : The frequency interval between two neighbouring points, which is $\Delta f = \frac{1}{T_0}$, is not used, nor computed by the fft. You should explicitly generate a frequency support vector if you want to plot the spectrum vs the right frequency for visualization purposes. Knowing $\Delta f$, the start and end value for the frequency axis in the bilateral representation are :

   $$f \in \left[ {\frac{{ - N}}{2}\Delta f,\frac{N}{2}\Delta f} \right] = \left[ {\frac{{ - f_s }}{2},\frac{{f_s }}{2}} \right]$$

A question (Food for thought)

In order to increase the spectral separation between two neigbouring spectral components, a physicist decides to double the sampling frequency $f_s$ of an original signal, while keeping the same total duration $T_0$. Explain why he is wrong by stating the consequences of $f_s$ doubling on :

1. the frequency resolution of the DFT result

2. the frequency range spanned by the DFT result.

Therefore, what is the right action he should have taken in order to increase the frequency resolution of the spectrum ?

Your answer here :

$\Delta f = \frac{f_s}{N}$ so if we double fs but keep the window T we also double N and therfore the frequency resoultion doesnt change.

However the freqency range is $\left[ {\frac{{ - f_s }}{2},\frac{{f_s }}{2}} \right]$ and is therefore doubled, but this wasnt intended.

In order to increase the freqnecy resolution the window T needs to be increased.

Exercise VII.2

1. Using the DFT (function fft , evaluate the Fourier transform of an even rectangular pulse $x(t) = A \ rect(\frac{t}{T_1})$ with a width $T_1=100 ms$ and an amplitude A=2, and plot it in the first graph of a multi-graph figure. Chose the right total sampling duration so that the frequency interval $\Delta f$ of the DFT equals 1 Hz.

   Use the instruction “stem” to plot the real and imaginary part of the DFT in the second and third graph of the multi-graph figure, and its modulus and phase in a third and fourth graph. The four graphs should be correctly indexed horizontally with the right frequencies, using a bilateral representation, i.e. showing both negative and positive frequencies. Don’t forget that in the raw DFT result, the negative frequencies are in the second half of the vector, and should thus be rearranged using e.g. fftshift.

2. As far as the amplitude of the modulus or the real part of the DFT at f=0 is concerned, do you find the right amplitude ? If not, didn’t you forget to multiply the DFT by the sampling period $\Delta t$ ?

For f = 0 we get $FT(0) = \frac{1}{N} \Sigma_n \, x_n$ we can compare this

3. What about the phase of the FT evaluated by the DFT ? Does it ensure that the FT is real, as should be the TF of a real and even signal ? Does this allow you to plot directly the FT (suppose it is contained in the vector X) in a single graph using a plot(X) or plot(f,X) ? If not, you probably obtain a plot of points in the complex plane; can you explain how the Hermitian property of the FS coefficients explains the observed symmetry of this graph ?

the imaginary part is very small and only non zero because of the computation. We can set these values to zero to ensure the DFT is real.

4. Verify that you have obtained a completely real FT by computing max(imag(X)), which should return exactly zero. In the negative, it might well be that your input sequence is not properly arranged. Remember that for the DFT, the input time sequence should start with the value taken by the signal at t=0, whereas the negative values are found in the second half of the input vector. In addition, if the first P points are at A=2, you should have only P-1 points at A=2 at the end of the input vector (hint : starting from a centered rectangular impulse, you can cleverly use the instruction fftshift in the time domain in order to arrange the input sequence properly ).

It seems tqht the condiotion is fulfilled but still we dont get exactly zero

5. Compare the effect of the following two modifications while keeping the same input signal :

   • Keeping the same total duration, double the number of points of the initial time sequence by doubling the initial sampling frequency. Describe the changes in the DFT output.

   • Keeping the initial sampling frequency, double the total length of the input sequence by introducing a number of zeros (a kind of “zero-padding” process) at an adequate place (Where is the right place ?), and describe what are the changes in the DFT output.

   • Which action is the most interesting as far as the frequency precision (resolution) of the DFT is concerned ? And what could be the usefulness of the other one ?

6. Use the DFT to evaluate the Fourier Transforms of two rectangular pulses with different durations $T_n$ (take $T_1=100 ms$ – it’s done already ! -, $T_2=50 ms$ and $T_3=200 ms$). Comment on the evolution of the width of the first lobe of the cardinal sines in the DFTs w.r.t. T~n~.

This is waht we expect increasing the width in time spice decreases the width in fourier space

VII.3 Evaluation of the spectrum of a realistic signal using the DFT - Windowing and “zero-padding” :

Background

The DFT may be used to evaluate the spectrum of a continuous-time signal x(t) after its analog-to-digital (ADC) conversion to x[k] at a suitable sampling frequency $f_s$. Such a signal is obviously a time-limited signal, and therefore we apply as above :

Exercise

The goal of this exercise is to use the DFT to estimate the spectrum of a « real » signal produced by a transducer. Load the file “experim.mat” into your environment. It contains a vector y with a discrete-time signal from a vibration transducer (also called « sonometer », although it is generally sensitive to spectral components outside the range of the human ear), sampled at a frequency of $1 kSa \ s^{-1}$ (kilo-sample per second).

1. Plot the graph of the signal y. Is it possible to recognize a priori some frequency components of this signal from the mere temporal representation ?

2. Use the DFT to compute the spectrum of this signal, and plot its modulus, taking care of the correct indexation of the (bilateral) frequency axis.

3. Windowing : Apply a Hanning window (function HANNING) to the temporal input signal before performing the spectrum evaluation with the DFT. Describe the main changes on the spectrum w.r.t. the spectrum of the original signal.

4. Zero-padding : The physicist suspects that the peak situated near 22 Hz could be a doublet made of two closely-spaced components. In the zero-padding trick, one adds a long sequence of zeros to the original input sequence, making it artificially longer. Do so by adding such a sequence so that the original duration will be multiplied by four, then evaluate the spectrum of this longer signal using the DFT. What happens to the peak near 22 Hz ? What are your conclusions ?

5. When you add such a zero sequence following some high values of the signal, there will be an abrupt jump of the long signal at the junction, which could appear on the spectrum as spurious high frequencies. Propose and apply another clever way to elongate the initial time sequence without adding any spurious frequencies ?

More details about the DFT

In this lab devoted to DFT implementation you can encounter some difficulties concerning the following points :

1. Obtaining purely real (or imaginary) $X_k$ coefficients (for cosine and sine power signals) in k=±1, therefore in $f=\pm f_c$.

2. The construction of the frequency scale.

3. Obtaining a flat phase when the module ||$X_k$||=0.

Obtaining purely imaginary or real $X_k$ for (co)sinusoids

The goal of the first exercise of this lab is to calculate the coefficients of the Fourier series development of a periodic signal, using digital, thus discrete-time signals. The condition of adequate windowing is to generate exactly one or m periods of a digital signal which is really periodic.

For a discrete-time sinusoid, the periodicity condition is that the discrete frequency $\nu_0$ in $x([n] = A \sin(2 \pi \nu_0 n))$ is a rational number $\nu_0 = \frac {k}{N} = \frac {1}{r}$ , and in this case the period is r points (in the figure below, r=5).

In our case, this is equivalent to imposing a condition on the sampling period, which must be a multiple of the frequency f~c~=1/T of the analog signal (see below).

So how do you get a periodic discrete-time signal ? First of all by avoiding to set the sampling frequency $f_s=1/T_e$ independently of the signal frequency $f_c$, as we can see below :

fc=3 ;
m=2 ;
te=0.031 ;	% not really a very good idea !
T0=m/fc ;
t=0 :te :T0-te ;
N=length(t) ;
X=1+2*sin(2*pi*fc*t) ;
…

It can be seen above that the digital signal is not periodic (note the shift of the points near the maximum), which gives rise to a non-zero real part in f=±3 Hz, which is associated with a phase, and to a frequency spread (non-zero modulus values other than in f=0 and f=±3 Hz) due to the discontinuity of the time connection.

In the example below, we set $f_c=3$, then $f_s=13 f_c$ from which we derive $t_s=\frac{1}{f_s}$ : you can see that the discrete sinusoid $x[n] = 1 + 2 \sin(2 \pi f_c n T_e) = 2 \sin(2 \pi \frac {n}{13})$ (from which $f_c$ has now disappeared, which ensures a frequency parameter $\nu = \frac{1}{13}$ , thus a rational ! ) is now rigorously periodic with a period n=13 (note that numerical values close to the maximum reproduce exactly from one period to another! ), and that $X_0$ is indeed the only coefficient with a non-zero real part.

fc=3
m=2
fe=13*fc;	%  A much better way !
te=1/fe;
T0=m/fc
t=0:te:T0-te;
N=length(t);
x=1+2*sin(2*pi*fc*t);
...

Construction of the frequency support

If we consider the raw result of the DFT X=fft(x) (without doing fftshift), the first of the N points of the vector corresponds to the frequency f=0, then we find the $X_k$ corresponding to the positive frequencies increasing with every multiple of $f=1/T_0$, the frequency step, up to half of the interval. Beyond that, we find the coefficients corresponding to the negative frequencies :

Thus, if N is odd, the coefficient $X_0$ is naturally found after the rearrangement due to the fftshift "in the middle" of the result vector, thus in X(1+N/2), with the same number of points on the right as on the left, and the range of the frequency scale is simply :

$f \in \left[ {\frac{{ - N}}{2}\Delta f,\frac{N}{2}\Delta f} \right] = \left[ {\frac{{ - f_s }}{2},\frac{{f_s }}{2}} \right]$

A small difficulty is introduced if N is even (much more favourable for the speed of the calculation of the fft !), since one cannot obtain the same number of points to the left and to the right of $X_0$. In this case, the result of the fftshift(fft(x)) is visible on the bottom right square of the following summary diagram (an example of DFT with N=10, assuming real $X_k$ :

We see that at the beginning of the vector fftshift(fft(x)), N/2 points correspond to negative frequencies, then $X_0$, followed by (N/2)-1 points which correspond to positive frequencies. There is thus one positive frequency less than the negative frequencies.

To reconstruct the right frequency scale regardless of the value of N (even or odd), we can use a parity test, which uses rem(x,y), "rem" for "remainder", which represents the remainder of the division of x by y. Thus, if rem(N,2) is null, it means that N is divisible by 2 without remainder, so N is a multiple of 2, so N is even. You can check that the instructions below allow you to find the right frequency scale:

if rem(N,2)~=0 % N odd
N2=(N-1)/2;
f=-N2:N2
f=f/T0;
else % N even
N2=N/2;
f=-N2:N2-1
f=f/T0;
end

Note : Avoid the linspace instruction which does not guarantee a zero crossing.

Obtaining a flat phase when the module |$X_k$|=0

When the module of a complex number is zero or almost zero, Octave/MatLab commonly displays any phases (although the odd character of the phases is preserved !) that may at first sight confuse the attentive observer :

However, when the modulus of a complex number $z = a + jb = \left\| z \right\|e^{j\theta }$ is zero, it can only be the origin of the complex plane, regardless of the phase of z (an analogy : when you are exactly at the North Pole of the earth, you can go around the earth, turning on yourself while staying there − $\theta$ merely varies between 0 and $2\pi$).

To get rid of these disturbing phase values, one solution is to make a test identifying the indexes of the coefficients of the TFD vector for which the modulus is zero, then to selectively set to zero the phases of the TFD vector for these indexes only, which gives a result more in line with the expected one :

Computing the spectrum of even signals

In computing the spectrum of even signals, it is especially important to make sure that the signal is well presented by 2N points (choosing $N=2^m$ increases the speed of calculation of the fft) to avoid any phase problem, remembering that :

• the first N points describe the positive times starting with t=0

• the last N points describe the negative times ending in -t~s~, with t~s~= sampling period make sure to generate a temporal signal which effectively passes through t=0

A simple possibility is to construct the time support vector of the signal as follows :

t=-N:N-1;
t=t*ts;
N=length(t);
x=5*sinc(t);        % or any other function generating 
					% the signal centered at t=0,
x=fftshift(x);

In short, the aim is to present the signal like this to the DFT, because only at this price will it interpret the signal as even :

Note that $x(1)=2$, while $x(end)=2-\delta$ ! The instruction X=te*fftshift(fft(x)); gives the following results, where we have a real and even DFT, in accordance with the general FT properties :

Conclusions

Complete this cell with your personal conclusions about what you've learned. Give an estimation on the time you've spent (outside the alloted time for lab sessions) realizing this lab session.