AOE3054 - Instrumentation Session 4
SPECTRAL ANALYSIS
A. Borgoltz, W. J. Devenport, W. L. Neu, S.R. Edwards, and N. Hari
Last revised
27 Mar 2022
In Experiment 6 analog, a function generator was manually controlled to excite a beam. The Instrumentation Week 5 Lab, Experiment 6 digital, is essentially a redo of the first Experiment 6, but will incorporate new digital measurement techniques to automate most of the data taking. Specifically, the myDAQs will be used to output a voltage signal that will control the function generator. The myDAQ will also measure the function generator output as well as the output from the proximeter. All operations will be controlled via Labview, using a code that builds off of the homework assignments and will be completed in Instrumentation Lab 4.
The objective of this lab period is to show you how to determine the frequency content of a signal measured using a computer. You will get to write data acquisition programs that include spectral analysis, and learn how to interpret spectra when you measure them, as well as avoid measurement pitfalls.
This chapter provides a fundamental background to the subject of spectral analysis including several LabView programs that you can run on your own computer to illustrate the techniques and problems. You will need to read and master this material before coming to lab. During the lab period you work in groups of two or three students to master this material and combine it with data acquisition under the direction and guidance of an instructor. You are expected to take notes (electronic or paper) to record your experiences. It is important that you save your notes as you will need them to prepare for and perform the digital dynamic response testing you will do on the experiment 6 rig during the final lab period.
You will be given
further specific instructions during lab. Copies of the lab material will be made available
on this page, linked at the bottom of this page and from the left hand column. The material below gives a
complete description of the techniques you will be using in lab. You will need to have read this manual before going through the lab slides.
Consider a signal generated by a probe (such as a hot-wire anemometer) sensing the velocity in the wake of a circular cylinder. The velocity fluctuates with time because the flow is unsteady.
 |
 |
Imagine you are measuring this signal as part of an experiment to understand the flow and its effects. The portion of the velocity signal shown clearly contains important frequency information. We see fluctuations at several frequencies apparently connected to the passage of eddies over the probe and the frequency at which eddies are shed from the cylinder. If the cylinder were supported by a flexible structure, knowing this shedding frequency would tell us the frequency at which that structure would be excited and thus help us analyze its structural response. So this information is clearly valuable for several reasons, but how do we get it? How do we determine what frequencies are in the signal?
When a signal contains only one frequency (i.e. it is sinusoidal) determining that frequency, its amplitude and phase can be a straightforward process (we discuss one such method at the end of this chapter). However, real signals are usually not sinusoidal, contain many frequencies and may contain random elements. Determining the frequency content of such a signal requires more sophisticated methods, referred to collectively as spectral analysis. The primary purpose of this chapter is to explain the methods of spectral analysis, their usage, capabilities and limitations.
Consider the general expression for a sinusoid, using the cosine function.
The function has a frequency f (in Hertz) that is equal to the inverse of the time it takes to complete one period of oscillation. At a given frequency it takes two pieces of information to specify such a wave; its amplitude A and its phase e (in radians). Note that the amplitude is half the peak to peak fluctuation and, if you work it out, the r.m.s. is 1/√2 (= 0.7071) of the amplitude. Note that the phase is an angle, and that negative phase corresponds to positive time delay of the wave.
When we ask the question "what frequencies are in the signal?" what we really mean is "If we fit the shape of the signal to a sum of sinusoids of different frequencies, what is the distribution of amplitudes and phases as a function of the frequency?" This question has a useful answer because of Fourier's theorem which states that any signal of zero mean value can be represented as a unique sum of sinusoids.
When we decompose our signal into its component sinusoids, the resulting plot representing amplitude or phase as a function of frequency is referred to as a spectrum. There are three types of spectra we need to be concerned with at this stage:
- The amplitude spectrum - a plot of the sine wave amplitude vs. frequency. Note that when engineers refer to the amplitude spectrum they may either mean the amplitude itself A or the r.m.s. amplitude (= 0.7071A). You often have to look carefully (or ask) as to which is being used.
- The phase spectrum - may be plotted in radians or degrees.
- The power spectrum - plot of Amplitude2/2 vs. Frequency. This is called the power spectrum because the square of the variable represented by the amplitude (e.g. velocity, voltage, displacement) is often proportional to a rate of work being done. However, the term power spectrum always means Amplitude2/2 vs. Frequency even when this isn't the case.
3. The Spectrum of a Sampled Signal
The idea of obtaining a spectrum from a measurement may seem overwhelming, not least because signals in the natural world can contain infinitely many frequencies. However, such continuous signals can also be broken into infinitely many time steps and we can measure their behavior in time by sampling them at regular intervals over some limited time. In an exactly analogous way, measuring a spectrum is an exercise in sampling it at regular intervals in frequency over a limited frequency range. To understand how this comes about we need to consider the whole measurement process.
Consider
once more the velocity signal from our cylinder wake. To measure this signal we
take samples of it at regular intervals Dt. We
already know that the maximum frequency we can observe in such a sampled signal
is the Nyquist frequency, equal to half the sampling rate. Thus by sampling in time
we already have set the highest frequency we can measure as 1/(2Dt).
(There may be problems, of course, if our original signal contains frequencies
higher than this, so our measurement is aliased, but more about that later.)
Inevitably we can only take samples of the signal for a finite time called the total sampling time T. We refer to this process of extracting a finite length piece of the signal as windowing (the idea being that we are looking at the signal through a window in time). The lowest frequency we can unambiguously infer from a windowed signal is one with a period that lasts the total sampling time. Thus, if N is the number of samples we took, then the total sampling time T = NDt and the lowest frequency we can measure is 1/T = 1/(NDt). Note that this also turns out to be the smallest difference in frequency we can infer from the windowed signal.
So when we sample the signal at time intervals of Dt over a total time T = NDt it turns out that we are also sampling its frequency content at intervals of 1/(NDt) up to a maximum frequency of 1/(2Dt). As a practical example, if our velocity signal is being sampled at a rate of 10 Hz (Dt = 0.1 s) and we are taking the 32 samples shown then we have sampled its frequency content at intervals of 0.31 Hz = 1/(32×0.1 s) up to a maximum frequency of 5 Hz=1/(2×0.1 s).
Now consider the problem of calculating the spectrum, i.e. expressing our sampled signal as a sum of sinusoids each of different amplitude and phase. Since the sinusoids will have a maximum frequency of 1/(2Dt) and a frequency interval of 1/(NDt), we will have 1/(2Dt) ÷ 1/(NDt) = N/2 of them. This makes good mathematical sense. Since we need two pieces of information (amplitude and phase) to define each sinusoid, we will have N pieces of information to obtain using the N samples of the signal that we originally took. Obtaining the mathematical expressions for the amplitude and phase in terms of the original sampled values of the signal can thus be thought of as an exercise in solving N equations for N unknowns. Performing the math, we find that the amplitude and phase of the nth sine wave, having a frequency fn= n/(NDt), is:
where the coefficients an and bn are given in terms of our original signal v(t) as:
Note that the index i in these expressions just references the samples that we took so, for example, i = 3 references the third sample we measured, v(3Dt). These expressions, for determining the amplitude and phase as a function of frequency, are referred to as the Discrete Fourier Transform. The discrete Fourier transform is closely related to the continuous Fourier transform used in analyzing linear systems and, for example, in controls and dynamic response problems. This close relationship makes the spectral analysis particularly useful when we use it in experiments that deal with such systems or problems.
The expression for our original sampled signal in terms of the sum of these sinusoids is:
This expression (which effectively reverses the process of determining the amplitude and phase from v(t)) is referred to as the Inverse Discrete Fourier Transform. Note that it contains one extra coefficient A0/2 that we didn't anticipate, apparently for zero frequency. Following through the expressions for an and bn for n = 0, however, you will find that this extra coefficient simply represents any mean value of the signal which cannot be accounted for in the sum of sinusoids.
The figure below shows spectra obtained by taking the discrete Fourier transform of our example velocity signal.
Note that the spectra are plotted against frequency index n, Interpreting these spectra is not hard if you remember the meaning of the frequency index n plotted on the horizontal axis. Specifically n = 1 corresponds to a wave with a wavelength that fits one time into the sampling window, n = 2 is a wave that fits twice into the sampling window, and so on. The power spectrum shows large amplitudes for waves that fit 3 and 6 times into the sampling window, with phases of about -40 and 45 degrees respectively. Look at the sampled and windowed signal and see if you can identify these components in the original data by eye. Confirm in your own mind that, for a 10 Hz sampling rate, these peaks are occurring at frequencies of 0.94 Hz and 1.88 Hz.
Finally, note that the phase values only have meaning when the amplitude is non zero (a sine wave with zero amplitude looks the same whatever the phase is). Thus, in the present case, the phase values for frequency indices of 11 to 16 means little.
In this, and the following section, we will look at some of the practical issues of computing spectra from a measured signal, with particular reference to MATLAB. We will examine a number of examples that involve spectral analysis of sine-wave signals. We use these single frequency signals because they are easily understood and therefore reveal clearly both the capabilities and limitations of spectral analysis. However, don't forget that the real power of spectral analysis is that it can be applied to any signal, whatever form it has, and however many frequencies it contains.
While the above equations for computing the discrete Fourier transform and its inverse are entirely correct, they are rarely used explicitly. This is because they contain a lot of redundant arithmetic. Much faster algorithms for computing the same results are usually used. These are called the Fast Fourier Transform (FFT) and Inverse Fast Fourier Transform (IFFT). We won't describe these algorithms here and you don't need to consider programming them since they have been part of standard software packages ever since such packages have existed (an indication of how important and universal spectral analysis is). Such packages often scale and present FFT results differently, however, so it is important that you find out what your software does before you use it. Below we look at computing FFTs in Matlab.
4.1 Matlab
Suppose the original samples of our signal in an array
v with elements 1
to N. The command c=fft(v);
computes the fast Fourier transform of v
and places it in the complex array c. To get the above coefficients from
c we have to scale it and separate it into real
and imaginary parts. We would write
a=real(2*i*c/N);
b=imag(2*i*c/N);
to get arrays a and b containing the coefficients an and bn (where the nth array element corresponds to the nth coefficient). Note that in Matlab i refers to the square root of -1. With the coefficients an and bn it is a straightforward matter to go ahead and calculate the amplitudes and phases and plot their spectra.
Performing the inverse transform basically means reversing this process. Starting with arrays a and b containing the coefficients an and bn we would write...
c=(a+i*b)*N/2/i;
v=ifft(c);
...to recover the sampled time signal v.
The ultimate objective of spectral analysis to determine the actual distribution of frequencies in a signal. This doesn't appear to be too difficult in the above examples, but these signals have been carefully chosen to avoid problems. Real signals may suffer from aliasing, may not fit well in the measurement window and may contain random elements or be corrupted by noise. Knowing how to handle these problems is key if you want your measured spectrum to be of some value.
5.1 Avoiding Aliasing Errors
In section 3 we discussed how the maximum frequency we
can observe in a sampled signal is the Nyquist frequency, equal to half the
sampling rate. What happens then if the original signal contains higher
frequencies? Do these frequencies simply not appear in the measured spectrum
since they are outside its range, or do they somehow corrupt the measured
spectrum causing error at the frequencies we can resolve? The second answer,
unfortunately, is the correct one.
One can write a very simple code to analyze the frequency content of a single sinusoid. In Matlab, such a code is reproduced below:
Let's assume we are looking at a 10 Hz sine wave being sampled at 1000 Hz for 1 seconds (1000 samples). To capture this signal accurately and avoid aliasing, we would need to sample at at least 20Hz (which is twice the frequency of the signal). At 1000Hz, we are well above that criterion and we see that the signal is captured accurately. Indeed the spectral analysis provides a measured frequency of 10Hz.
But if we sample the input acquisition signals at a much lower rate, we will see that the number of points representing the signal reduces and the signal becomes more and more inaccurate.
Indeed, for a lower input sampling rate of 30 Hz (and 30 samples to maintain the 1s smapling period), we can see that the signals captured becomes inaccurate. Indeed, this looks more like a triangular wave than a sine wave. We are now just above the Nyquist requirement of 20Hz, which means the signal is not aliased. The figure below shows that the number of points in the acquired signals is too low to provide a faithful estimate of the original signal. However, the spectral analysis does provide an accurate estimate of the frequency content.
If we use an even lower sampling rate of 15Hz which is much smaller than the Nyquist criterion here, both the amplitude and the frequency content of the signal are incorrectly captured. Indeed the frequency estimation using FFT gives an inaccurate value (5Hz compared to the signal frequency of 10Hz). This is because the sampled data is aliased and the frequency of the original sine wave cannot be accurately captured.
In general aliasing creates two problems. First, it can cause us to mistake a high frequency signal for a low frequency one. This is common, for example, when electrical interference from a computer contaminates a sensor signal. The interference (usually a signal at many MHz) can get aliased and appear as though it is part of the legitimate low frequency sensor signal, and therefore be misinterpreted as a fluctuation in velocity, temperature, or whatever physical property we are trying to measure. Second, when a signal contains a broad range of frequencies, aliasing can corrupt the entire shape of the spectrum.
So how can we minimize the influence of aliasing on a spectral measurements? There are three measures we can take.
- Take samples at a higher rate, so the Nyquist frequency exceeds the highest frequency in the signal. This is ideal, but not always possible since all A/Ds have maximum sampling rates. Also you may not know what the highest frequency is. Signs that your sampling rate may be high enough are (i) that the spectral values you calculate at the highest frequencies are zero or almost zero (only works for broadband signals), or (ii) that spectra of the same signal measured at different sampling rates have the same shape.
- Use an analog electrical filter to remove all the frequencies in the signal above the Nyquist before you sample it. (While you don't get the whole spectrum in this case, at least the part of the spectrum you do get is uncorrupted).
- Be aware of the problem and alert to its possible presence. You may have to accept some minimal level of aliasing in your measurement. If you are aware that it is there, then you won't mis-interpret its effects.
5.2 Spectral Resolution
Broadening of the spectrum describes what happens when a signal does not fit neatly into the measurement window.
The figure below shows the output from the spectral analysis code written during the lecture (and during Session 4)
in its default configuration (that shows an amplitude spectrum computed from 100 samples, taken at a rate of 1 kHz,
of a 50 Hz sine wave with an amplitude of 1 and zero phase).
The spectrum looks perfect, but this is partly due to a lucky choice of parameters. The numerical Fourier transform treats the measured part of the signal as though it were part of an infinitely extending signal made up of repetitions of the measured part. At 50 Hz exactly 5 periods of the wave fit into the measurement window between t = 0 and 0.1 s, so this measured part of the signal exactly fits to its copy from 0.1 to 0.2 s which, in turn, fits to the copy from 0.2 to 0.3 s, and so on. The implied infinitely extended signal is just a continuation of the 50 Hz wave we measured, resulting in the perfect spectrum.
Of course this doesn't often happen in real life. The chances of picking a total measurement time that is an exact number of periods of the oscillation of a structure, or the flow rate in a fuel system are almost zero. Indeed, many signals we encounter in engineering measurement won't be periodic at all. So what happens in this case? Let's re-run the spectral analysis code but with a signal frequency to 45 Hz.
Now there is an extra half period at the end of the measurement window, so the implied repeated signal has a jump at every repetition as illustrated below:
This rather ugly signal implies a range of frequencies (not just 45 Hz) and the spectrum we get reflects exactly that. There is still a peak around 45 Hz, but its amplitude is reduced and there are significant spectral levels over a whole range of surrounding frequencies. This is the problem of broadening.
When specifying your samping scheme, it is therefore important to understand the impact the scheme has on the information that can be extarcted from the samples signal. When selecting a sampling scheme, one should make sure to note the following:
- We know that the sampling rate of the input channel should be higher than double the frequency of the excitation signal (Nyquist criterion). Remember that this is a minimum. Typically you will want a much higher value (1000 Hz is a good start for the beam setup) to capture more data points and define the signal better.
- If we want fine increments in the measured frequency (maybe to narrow down the estimates of the natural or resonance frequency), then we need to make sure that the sampling period is long enough to capture these small increments. If we want 0.5Hz, we would need 2s of acquisition, 0.1Hz increments imply 10s etc...
For instance, let's look at the figure below, where three cases are presented. For each case, the resolution of the spectrum df(the smallest measurable change in frequency) implied by the sampling scheme is provided:
- Case 1: 1000 samples, taken at a rate of 1 kHz, of a 13 Hz sine wave giving a df = 1 Hz
- Case 2: 500 samples, taken at a rate of 1 kHz for the same 13 Hz sine wave gives a df = 2 Hz
- Case 3: 100 samples, taken at a rate of 1 kHz, gives a df = 10 Hz.
Essentially, the frequency resolution (df, smallest increment in the spectra) must be the inverse of the sampling period:
Note that broadening doesn't affect signals that fall to zero at the ends of the measurement window (such as measurements from short duration phenomena like transients). Nor does it have much impact on signals where the amplitude spectrum varies smoothly from frequency to frequency, and there are no sharp spikes. However, when we have a signal containing one or more dominant frequencies (as in the present example) broadening can result in significant smoothing of the sharp spikes that would otherwise appear in the spectrum.
In this case the effects of broadening can be minimized by multiplying the signal by a smooth function that goes to zero or nearly zero at the ends of the record, called a window function. There are a variety of different commonly used window functions. Some are listed in the table below and illustrated in the following figure. All these have the form
where T is the total record length and the constants c0 , c1 , and c2 are as follows:
| Window type | c0 | c1 | c2 |
| Blackman | 1.0000 | 1.1905 | 0.1905 |
| Blackman-Harris | 0.9790 | 1.1509 | 0.1833 |
| Exact Blackman | 1.0000 | 1.1640 | 0.1801 |
| Hamming | 1.0000 | 0.8519 | 0.0000 |
| Hanning | 1.0000 | 1.0000 | 0.0000 |
5.3 Dealing with Randomness
There are two kinds of randomness we have to deal with in measurements;
unwanted noise that masks a sensor signal, and thus the variations in the
physical quantity we are trying to detect, or randomness in the physical
quantity itself, such as in the velocity of a turbulent flow or in the unsteady
lift of a stalling wing. Both of these situations call for averaging, but of different spectral quantities.
We can eliminate unwanted noise from a repeating signal by measuring multiple repetitions, and then averaging them together before taking the spectrum. As long as each repetition of the signal appears at the same position in the window each time we measure it, then the repeated signal will average to its true value, whereas the random noise fluctuations will average to zero. The same result can be achieved by calculating the spectrum of each repetition, and then averaging the raw spectral coefficients an and bn (before calculating amplitude and phase).
This exact situation is encountered when using computer data acquisition to measure the dynamic response of a structure. The response can be revealed for all frequencies, for example, by measuring the oscillations in time of the structure after it experiences an impulsive force (e.g. it is tapped with a hammer) and then calculating the spectrum (more about this later). There will be noise in the sensor not just from interference and electronics but also in the movement the structure picks up from other sources, such as building vibrations. So if the response is determined from a single measurement it will have uncertainties associated with this noise. Instead, the way to proceed is to repeat the measurement multiple times, triggering the data acquisition to start at the same instant (when the hammer hits, perhaps) each time, and then averaging the spectral coefficients. The more repetitions (averages) the lower the uncertainty in the spectrum will be.
Sometimes the random part of a signal is just as important as any deterministic part. The flow velocity in the wake of a circular cylinder high speed may often consist of a periodic fluctuation associated with vortex shedding, super-imposed on random turbulence. In this circumstance we will likely want to know the spectrum of the turbulence as much as we want to know that of the vortex shedding. Another example can be found in measuring the dynamic response of a structure. A second way to get the response at all frequencies is to measure the spectrum of the random movements produced by the structure in response to randomly fluctuating force (this technique will be discussed in more detail later).
To obtain the average spectrum of a signal with random components, we average the amplitude squared (An2) at each frequency, instead of the Fourier coefficients. This way the negative and positive fluctuations in the amplitude, produced by the randomness, don't cancel out. When the averaging is complete we can estimate the amplitude spectrum by taking the square root of these averages and the power spectrum as half these averages. Note that this type of averaging eliminates any phase information.
During your instrumentation session, you will go through a set of slides (to be posted on CANVAS) that will guide you on how to implement advanced programming options in an automated LabVIEW code for Experiment 6.