Essential MATLAB for Engineers and Scientists (2013)
PART 2
Applications
CHAPTER 15
Signal Processing
Abstract
The objective of this chapter is to introduce signal processing. We apply two methods, of the many available methods, to examine large data sets. They are harmonic analysis and fast-Fourier transform (FFT) analysis.
Keywords
Sines and cosines series; Fourier expansion; Harmonic analysis; FFT
CHAPTER OUTLINE
Harmonic analysis
Fast Fourier Transform (FFT)
THE OBJECTIVE OF THIS CHAPTER IS TO INTRODUCE SIGNAL PROCESSING. YOU WILL LEARN ABOUT TWO METHODS, OF MANY, TO EXAMINE LARGE DATA SETS. THEY ARE:
■ Harmonic analysis.
■ Fast Fourier Transform analysis.
In this chapter we apply MATLAB to examine periodic data sets. Harmonic analysis is introduced in the next section to examine periodic data sets. In the subsequent section the fast Fourier transform (FFT) method is applied to examine the temporal-evolution of measured or computed data. These procedures are important methods in the field of signal processing.
A signal is a sequence of data that is measured by experimental methods or computed by computational methods. For example, the data 
is measured at relatively small intervals of a particular variable x, e.g., x could be time or a space dimension. If the signal is periodic in a range of x, e.g., −π to π, harmonic analysis is a useful and quite practical technique. The FFT method is a related method. Both methods are associated with the fitting of data with a Fourier series. In addition, both methods are widely used in engineering, e.g., in naval architecture, acoustics (or noise) and vibration engineering, electrical engineering, ocean engineering, etc.
The theory of Fourier series is an important topic in mathematics for physicists and engineers. A reasonably thorough discussion of this topic can be found in Mathematics of Physics and Modern Engineering, Second Edition, by I.S. Sokolnikoff and R.M. Redheffer, a book published by McGraw-Hill (1966). Because of the importance of this topic, it should not be surprising that there are many other books that describe the theory behind the method and the practical applications of the method. In this introductory section we give the results associated with applying the Fourier series method that are needed to examine signals.
Let us suppose is given by the trigonometric series
(15.1)
where 
and, hence,
(15.2)
for 
, 1, 2, …, and
(15.3)
for 
, 2, 3, …. This assumes that the function is periodic, i.e., 
.
If we are given a data set and it is assumed to be periodic in the interval , we can compute the coefficients 
and 
by applying the above formulas. With this approach we say that we expanded in or say that has the Fourier series
(15.4)
where N is a positive integer greater than or equal to 1. It is the number of terms to be evaluated to construct a Fourier series (or Fourier expansion) of . If 
and is continuous and periodic, the last equation is essentially the trigonometric expansion of and the result is exact. An important advantage of this series is that it can represent discontinuous functions as well. Of course, in this case or in cases where the trigonometric representation is not exact, it is a best fit representation of in terms of the assumed periodic function, Equation (15.4).
15.1. Harmonic analysis
Let us examine the data illustrated in Figure 15.1 by fitting a trigonometric series through the data. This sample set of “observed” data was generated by using the following formula:
(15.5)
where the error was generated by the random number generator randn, a built-in function in MATLAB. Equation (15.5) was selected for this example because it is a five-term trigonometric serious plus an error term. It was selected so that we can examine how well harmonic analysis plucks out the harmonics (i.e., the first five terms) from a signal that contains error (the error in some applications of signal processing is also called “noise”).
FIGURE 15.1 Sample of ten repeated periodic data sets with error.
We want to fit a finite trigonometric sum to a set of observed data. In this case the discrete set of observed data contains the values of at 
, 
, 
, …, 
. It is assumed that 
. Note that 
is the domain of x that the data are observed. The values of are evaluated (or observed) at equally spaced intervals of 
. The example numerical analysis is for 
where 
is the number of equally spaced values of f from 0 to π where the data are measured. The data obtained at these points were computed by applying Equation (15.5). Since f is prescribed on the interval and it is assumed that , harmonic analysis is a good way to examine the data to determine its harmonic content, i.e., to determine the coefficients of the Fourier series that we are applying to expand the data. Note that the noise in the data introduced by the random number generator introduces an error in the periodicity assumption. Hence, in the example below ten sets of the observed data are examined by repeating the computational experiment ten times. The results illustrate that the error (or noise), in this case, is random and can be averaged away if a sufficient number of the same experiment is performed.
It turns out that if you have 
equally spaced data points we can fit the data with 
terms of a Fourier series. The results are exact in the sense that we have a system of linear equations in unknowns. It is reasonably well known (see, e.g., the book by Sokolnikoff and Redheffer cited above) that the solution to this problem is
(15.6)
(Note that 
is twice the value of 
in the Fourier series given in the introduction to this chapter.) In this case, it can be shown that
(15.7)
(15.8)
(15.9)
(15.10)
This system of equations was applied to examine the data generated by Equation (15.5). The MATLAB script developed and applied is as follows:
% Harmonic analysis
clear;clc
for ir = 1:10
N = 200; % Nc = N+1 implies x = 0 (center of domain).
% Np = 2N+1 are the number of points on the
% closed interval x = [-pi, pi].
dx = pi/N; % Spacing between the points on x.
x = -pi:dx:pi;
f = [.6 - cos(x) - .5*cos(2*x) + .4*cos(3*x) ...
+ sin(x) + 0.1*randn(1,length(x))];
A0 = (1/2/N) * sum(f(1:end-1));
A(N) = (1/2/N) * sum(f(1:end-1).*cos(N*x(1:end-1)));
for jn = 1:N-1
A(jn) = (1/N)*sum(f(1:end-1).*cos(jn*x(1:end-1)));
B(jn) = (1/N)*sum(f(1:end-1).*sin(jn*x(1:end-1)));
end
AA = [A0 A];
figure(1)
hold on
stem([0:length(AA(1:11))-1],AA(1:11),'ko','LineWidth',2)
stem(B(1:9),'-.kd','LineWidth',2)
legend('A_n''s', 'B_n''s')
xlabel('Mode number n')
ylabel('A_n''s or B_n''s')
title('Harmonic analysis example')
hold off
figure(2)
hold on
plot(x,f,'k')
title('Ten sets of data with error compared')
xlabel('x'),ylabel('f')
hold off
end
Note that the algorithm was executed ten times to examine the influence of the random error on a converged solution. The random error is plus-or-minus 10% of unity. The ten samples of the periodic data examined in the script is illustrated in Figure 15.1. The ten associated harmonic-analysis results are given in Figure 15.2. As we might have expected, the results, on average, are pretty close to the exact results described next.
FIGURE 15.2 Harmonic analysis of the ten periodic data sets, f, compared.
Since the function without error is
it is a five term trigonometric series. The harmonic analysis of this function is exact. This function and the exact results of the harmonic analysis of it are illustrated in Figure 15.3 and Figure 15.4, respectively. The figures were generated with the same MATLAB script given above with the following changes. The coefficient of the error term was set to zero as follows: We changed 0.1*randn(1,length(x))]; to 0*randn(1,length(x))];. Since the ten repeated trials are identical, if you change the repeated-trials for loop, for ir = 1:10, to for ir = 1:1, the results you obtain will be identical to the results in the figures for the zero error case.
FIGURE 15.3 Sample of the periodic data set, f(x), without error.
FIGURE 15.4 Harmonic analysis of the periodic data, f(x), without error.
Comparing the error-less data with the noisy data generated by adding a random error, illustrates the need for making repeated measurements in the laboratory to get a sensible result for the amplitudes of the harmonics that represent the shape of the phenomenon sought, viz., in the interval . Note that if a trigonometric series is used to define the input data set, the output of harmonic analysis yields an exact reproduction of the input series.
15.2. Fast Fourier Transform (FFT)
James W. Cooley and John W. Tukey published An Algorithm for the Machine Calculation of Complex Fourier Series in 1965; it appeared in Mathematics of Computations 19, 297–301. This paper was cited in the book edited by Leo L. Beranek and István L. Vér titled Noise and Vibration Control Engineering published by John Wiley & Sons (1992). This book provides a good review of signal processing as well as its practical applications in noise and vibration engineering. The discussion to follow follows closely the discussion in this book.
The Fourier transform of a time history signal, 
, that is measured over the time interval 
is defined for all frequencies, both positive and negative, by
(15.11)
where Y is a function of frequency, ω. In terms of a digital time series of N data values where 
, , 1, 2, …, 
. The Fourier transform may be written as a complex Fourier series
(15.12)
where the spectral components, i.e., the 
s, are generally complex valued and are defined only at N discrete frequencies
Hence, this is essentially the complex version of the Fourier series expansion introduced in the introduction of this chapter. To clarify this point further, let us divide Equation (15.12) by 
. This yields the following result:
(15.13)
Note the analogy between this formula and the formulas for the Fourier coefficients given in the previous section on harmonic analysis. In this case, however, the coefficients are complex numbers. Since MATLAB deals with complex numbers and does complex arithmetic, applying the FFT methodology to examine signals is reasonably straightforward. This will be illustrated below by an example.
The complex Fourier coefficients defined by Equation (15.13) are for an assumed periodic function and, hence, the time history is assumed to repeat over periods equal to the sampling period, T. The Fourier components are unique only out to 
, i.e., out to the frequency 
, commonly known as the Nyquist frequency, 
, of the digital signal. At this frequency there are only two sample values per cycle and, hence, the error known as aliasing is initiated; for details see, e.g., Beranek and Vér (1992) cited above. The first 
Fourier coefficients, from 
to , define the spectral components (i.e., magnitudes of the Fourier coefficients of significant magnitude) at positive frequencies, while the last 
Fourier coefficients, from 
to 
, define the spectral components at negative frequencies.
The various FFT algorithms developed since the publication of Cooley-Tukey algorithm are well documented in the technical literature. The details of developing and FFT algorithm is also beyond the scope of this presentation. However, since the procedure applied in the example examined in this section can be applied to other problems that may be of interest to the reader, a few of the hints that need to be considered in applying the Cooley-Tukey algorithm pointed out by Beranek and Vér are worth noting since this algorithm is the basis for the FFT algorithm applied in the example discussed in this section. They are:
■ It is useful to restrict the number of equally spaced data points for each FFT to a power of 2, i.e., 
, where 
to 12 are commonly used.
■ The frequency resolution of the Fourier components will be 
.
■ The Nyquist frequency where aliasing is initiated occurs at and it is 
. The sampling frequency should be at least twice the highest frequency expected to be found in the signal.
■ The first Fourier components up to the Nyquist frequency are related to the last components above the Nyquist frequency by 
, , 1, 2, … , , where the asterisk denotes complex conjugate.
■ The Fourier components that define only the positive frequencies is called a one-sided spectrum. This spectrum is given by 
, 
, and 
, 
, 2, … , .
The one-sided spectral components for a periodic signal with zero mean are given by
(15.14)
where 
is defined in Equation (15.12). The magnitudes of the Fourier components, 
, are usually plotted as a stem plot (see example below) and, hence, it is often called a line spectrum. Because each of the Fourier components are complex, they define both a phase and an amplitude for each of the components. However, the phase information is generally utilized only in applications where there may be a need to reconstruct the time history of the signal or to determine the peak values of the signal.
Let us consider the following example to illustrate how MATLAB can be used to apply an FFT to a signal. Consider the following signal:
(15.15)
where the sample is taken over the interval of time . The sampling frequency is 
. Thus 
. The power spectrum is also usually plotted to isolate the dominant frequencies in a signal. It is 
, where 
and 
was selected in this example. With this brief description of a signal processing problem, the following script was applied that implements the FFT methodology to examine this signal (it is based on an example in the help documentation available within MATLAB).
% ---------Example of the application of FFT------------
% Example of the application of FFT
Fs = 3000; % Sampling frequency
T = 1/Fs; % Sample time
pwr2 = 11;
L = power(2,pwr2); % Length of signal
t = (0:L-1)*T; % Time vector for FFT
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t)); % Sinusoids plus noise at t
plot(Fs*t/1000,y,'k'),axis([0 2 -10 8])
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (seconds)'),ylabel('signal y(t)')
figure(2)
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
Pyy = Y.*conj(Y)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
%
% ---------Plot single-sided amplitude spectrum.--------
FF = 2*abs(Y(1:NFFT/2+1));
stem(f(1,1:110),FF(1,1:110),'k','Linewidth',2)
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(\omega)|')
figure(3)
stem(f(1:110),Pyy(1:110),'k','Linewidth',2)
title('Power spectral density')
xlabel('Frequency (Hz)')
ylabel('Pyy')
The signal to be analyzed by spectral analysis (as it is called) is illustrated in Figure 15.5. The line spectrum and the power spectrum predicted for this signal are given in Figure 15.6 and Figure 15.7, respectively. There are two things to note. One is that the dominant frequencies in the signal are identified clearly in the spectra. The second is that the power spectrum isolates the dominant modes in the signal more clearly.
FIGURE 15.5 Example of signal to be analyzed by spectral analysis.
FIGURE 15.6 Line spectrum of the signal in Figure 15.5.
FIGURE 15.7 Power spectrum of the signal in Figure 15.5.