- Email: [email protected]
Frequency Analysis Part I
C H A P T E R
11 Frequency Analysis Part I: Fourier Decomposition
11.1 GOALS OF THIS CHAPTER This chapter introduces the most common method of decomposing a time series into frequency components, Fourier analysis. You will learn about the Fourier transform and the associated amplitude and phase spectra. The MATLAB® implementation of the fast Fourier transform (FFT), an efficient algorithm for calculating Fourier transformations, will be introduced and applied to the analysis of human speech sounds.
11.2 BACKGROUND Figure 11.1 shows typical recordings of two human vowel sounds. How can you characterize these different sounds? Frequency analysis provides a way to examine the relative contributions of various frequencies to an overall signal. In the case of an auditory signal, a given frequency component would be termed pitch.
11.2.1 Real Fourier Series Take some continuous function f. We can approximate such a function with a weighted series of sinusoids of various frequencies. Such a series is termed the real Fourier series: fðtÞ 5
N N X X a0 1 an cosðntÞ 1 bn sinðntÞ 2 n51 n51
ð11:1Þ
Here, the coefficients an and bn represent the relative strength of each frequency component n/2π. [a0 represents the nonoscillatory component of f(t).] So, given f(t), determining the coefficients an and bn allows for the representation of f(t) as a series sum of sinusoids.
MATLAB® for Neuroscientists. DOI: http://dx.doi.org/10.1016/B978-0-12-383836-0.00011-4
229
© 2014 Elsevier Inc. All rights reserved.
230
11. FREQUENCY ANALYSIS PART I: FOURIER DECOMPOSITION
0.2 0.15 Amplitude
0.1 0.05 0 −0.05 −0.1 −0.15
0
0.5
1
1.5
2
2.5
3
3.5
2.5
3
3.5
Time, in secs 0.15 0.1 Amplitude
0.05 0 −0.05 −0.1 −0.15 −0.2
0
0.5
1
1.5
2
Time, in secs
FIGURE 11.1
Acoustic time series representing two different human vowel sounds.
We will exploit two special properties of the sine and cosine functions to find the Fourier series coefficients an and bn. Over the interval 2 π to π, cosine and sine functions with differing frequencies have the special property of orthogonality. The integral of the product of two mutually orthogonal functions evaluates to zero. So, the integral of the product of cosine or sine functions with differing frequencies results in zero over this interval. Another interesting property of sine and cosine is that the integral of the square of a cosine or sine function over this integral is π. To find the strength, am, of a cosine component m, multiply by the corresponding cosine function and integrate: ðπ
ðπ fðtÞ cos ðmtÞdt5
2π
2π
N ð N ð X X a0 cos ðmtÞdt1 cos ðmtÞan cos ðntÞdt1 cos ðmtÞbn sin ðntÞdt 2 n51 n51 π
π
2π
2π
ð11:2Þ All terms on the right side except the cosine term where m 5 n yield zero: ðπ
ðπ fðtÞ cos ðmtÞdt 5 am 2π
cos2 ðmtÞdt 2π
III. DATA ANALYSIS WITH MATLAB
ð11:3Þ
231
11.3 EXERCISES
The right side integral evaluates to one over the integration range, yielding an expression for the Fourier series term coefficient: ðπ fðtÞ cos ðmtÞdt 5 πam
ð11:4Þ
2π
am 5
1 π
ðπ fðtÞ cos ðmtÞdt
ð11:5Þ
2π
In general, the interval of f(t) will not be 2 π to π. For an interval centered on x with length 2L, the expression becomes 1 am 5 L
x1L ð
fðtÞ cos
π L
mt dt
ð11:6Þ
x2L
A similar procedure using sine functions yields the coefficients for the sine terms of the Fourier series.
11.3 EXERCISES
EXERCISE 11.1 Write a MATLAB function to calculate coefficients for a real Fourier transform. Hint: The function will need to shift the
interval so that the interval encompasses the entire time series. In other words, x 5 0 and L 5 half the range of t.
11.3.1 Complex Fourier Transform Euler’s identity, eiωt 5 cos ωt 1 i sin ωt
ð11:7Þ
provides a straightforward way to formulate complex Fourier series representation for a given function, f(t): fðtÞ 5
N X
cn eint
ð11:8Þ
n52N
Similar to the real transform, coefficients for the complex Fourier transform can be found by
III. DATA ANALYSIS WITH MATLAB
232
11. FREQUENCY ANALYSIS PART I: FOURIER DECOMPOSITION
1 cm 5 2π
ðπ
fðtÞe2imt dt
ð11:9Þ
2π
for a given coefficient m over the interval 2 π to π. Over the interval x 2 L to x 1 L, this becomes 1 Cm 5 2L
x1L ð
fðtÞe2ðiπmt=L Þdt
ð11:10Þ
x2L
EXERCISE 11.2 Write a MATLAB function to calculate coefficients for a complex Fourier transform. This is essentially the discrete Fourier transform (DFT): Fk 5
X 1 N21 2πn fn e2i N k N n50
where k ranges from 0 to N 2 1, N is the number of points, and fn is the value of the function at point n.
ð11:11Þ
Let’s look at how this method of the Fourier transform scales with N. Given a time series with N values, this method requires a multiplication of the series and the corresponding Fourier component and subsequent sum for each coefficient. Assuming a number of coefficients equivalent to N, then you have a process that scales with N2. In other words, as N increases, the time required to compute the Fourier transform increases as N2.
11.3.2 Fast Fourier Transform With a few special tricks, a faster algorithm, the fast Fourier transform (FFT) that scales in N log N time can be formulated. One of these tricks involves taking advantage of datasets exactly 2N elements long. The increase in processing speed has made the FFT ubiquitous in signal processing. While a complete derivation of the algorithm is beyond the scope of this book, invoking the MATLAB implementation of the FFT will be discussed. MATLAB provides an FFT function fft(X), where X is a vector in time space. fft returns the frequency space representation of X. To visualize the importance of the difference in scaling, execute the following code: figure hold on N 5 1:10 * 100; plot(N, N.^ 2, 'b') plot(N, N.*log(N), 'r')
III. DATA ANALYSIS WITH MATLAB
11.3 EXERCISES
233
EXERCISE 7.3 If N represents sample size, what can you observe about the benefits of scaling as N grows? Where does the efficiency of the
FFT algorithm benefit most, for large N or small N?
11.3.3 The Inverse DFT As you might imagine, there is an inverse to the DFT: fn 5
N 21 X
2πn
Fk ei N k
ð11:12Þ
k50
MATLAB provides ifft() to perform the inverse discrete Fourier transform. EXERCISE 7.4 Generate a single sine wave. Use fft() to generate the discrete Fourier transform. Use
ifft() to retrieve the original sine wave from the DFT.
11.3.4 Amplitude Spectrum Often when you are using Fourier analysis, the amplitude spectrum is one of the first analyses performed. The amplitude spectrum graphs amplitude against frequency. In terms of the Fourier series representation, the amplitude spectrum depicts the magnitude of the coefficients at various frequencies. As such, it depicts the relative strengths of the various frequency components. The following code generates a time series composed of 10 sine waves whose frequencies and amplitudes vary systematically. L 5 1000; X 5 zeros(1,L); sampling_interval 5 0.1; t 5 (1:L) * sampling_interval; for N 5 1:10 X 5 X 1 N * sin (N*pi*t); end plot(t, X); Y 5 fft(X)/L; Now, the variable Y contains the normalized FFT of X. Note the normalization factor L. Displaying the amplitude spectrum of X requires plotting the amplitudes at various
III. DATA ANALYSIS WITH MATLAB
234
11. FREQUENCY ANALYSIS PART I: FOURIER DECOMPOSITION
frequencies. Note that fft returns only a single value, the transform coefficients. Now, how do you determine the frequency scale? The return value of the FFT assumes that frequency is evenly spaced, from 0 to a theoretical result called the Nyquist limit. Nyquist demonstrated that a discrete sampling of a continuous process can capture frequencies no higher than half the sampling frequency. Since the code above has the sampling interval, this Nyquist limit is half the inverse of the sampling interval. The following code calculates the Nyquist limit for the time series: NyLimit 5 (1 / sampling_interval)/ 2; When viewing the FFT, it is important to remember that the result is the complex transform. Thus, simply using the result of the FFT as a set of real coefficients can cause a number of problems. To display the amplitude spectrum, the absolute value of the complex coefficients will be used. The values returned by fft are the coefficients for frequencies from the negative Nyquist limit to the positive Nyquist limit. If the time series data are purely real, then the resultant transform will have even symmetry. That is, the transform will be symmetrical across the abscissa. So, in this very frequent case, only the first half of the result of fft is used. The following code employs linspace to generate frequency values and plots the amplitude spectrum. linspace generates a linearly spaced sequence of values given initial and final values. Here, the initial and final values are 0 and 1, with a value count of L/2. The resultant vector is scaled by the Nyquist limit to generate the frequency vector. F 5 linspace(0,1,L/2)*NyLimit; plot(F, abs(Y(1:L/2)));
11.3.5 Power Power at a given frequency is defined as ΦðωÞ 5 jFðωÞj2 5 FðωÞF ðωÞ
ð11:13Þ
where F* is the complex conjugate of F. To do this in MATLAB, use the function conj to return the complex conjugate of a series of complex values. Here is a plot of the power spectrum of the time series generated for the amplitude spectrum: plot(F, (Y(1:L/2).*conj(Y(1:L/2)));
11.3.6 Phase Analysis and Coherence A power spectrum alone is not a complete representation of the information in the original signal. The various Fourier components can have various phases relative to one another, as illustrated in Figure 11.2. You can plot relative phase by frequency by plotting the inverse tangent of the ratio between the imaginary component and the real component. Why is this the case? Imagine
III. DATA ANALYSIS WITH MATLAB
235
11.3 EXERCISES
FIGURE 11.2
Phase difference between two sinusoids with the same frequency.
2 sin(x) sin(x+1)
1.5 1 0.5 0 −0.5 −1 −1.5 −2
0
1
2
3
4
5
6
7
8
9
10
the complex plane, with pure real values along the abscissa (x-axis) and pure imaginary values along the ordinate (y-axis). Any complex value in your 1D Fourier transform can be represented with a coordinate pair. The magnitude of the value is simply the distance from the origin to the coordinates, or the complex modulus. The phase is the angle formed by the abscissa and the line passing through theorigin and the complex point. Thus, using basic trigonometry, the phase angle is tan21 imag real . How can you represent this in MATLAB? L 5 1000; X 5 zeros(1,L); sampling_interval 5 0.1; t 5 (1:L) * sampling_interval; for N 5 1:10 X 5 X 1 N * sin (N*pi*t); end plot(t, X); Y 5 fft(X)/L; phi 5 atan(imag(Y)./real(Y)); F 5 linspace(0,1,L/2)*NyLimit; plot(F, phi(1:L/2)); EXERCISE 11.5 Compare the phase spectrum generated in the preceding exercises with the phase
spectrum of the corresponding cosine function. Compare their power spectra.
III. DATA ANALYSIS WITH MATLAB
236
11. FREQUENCY ANALYSIS PART I: FOURIER DECOMPOSITION
TABLE 11.1 Average First and Second Formant Frequencies for Selected American English Vowels Vowel sound
First formant
Second formant
Bit
342
2322
But
623
1200
Bat
588
1952
Boot
378
997
(Data from Hillenbrand et al., 1995.)
11.4 PROJECT In this project, you will be asked to use Fourier decomposition to analyze vowel sounds produced by human speakers. On the companion website, you will find five examples of vowel sounds as produced by male American English speakers. Each sound corresponds to one of the vowel sounds in Table 11.1. The formant frequencies in Table 11.1 note the average formant frequencies as spoken by a male speaker of American English. You will use power spectra of these sounds to classify the recordings as one of these vowel sounds in the table. To complete this project, you need to understand how formants relate to frequency analysis. The human vocal tract has multiple cavities in which speech sounds resonate. As such, most sounds have multiple strong frequency components. In classifying speech sounds, the lowest strong frequency band is termed the first formant. The next highest is termed the second formant, and so on. Vowels lend themselves to a particularly simple characterization through their formants. Typically, vowel sounds have distinguishable first and second formants. Table 11.1 shows first and second formants for four vowel sounds in American English. Thus, the short “i” sound would have strong frequency representation at 342 Hz and at 2322 Hz.
MATLAB FUNCTIONS, COMMANDS, AND OPERATORS COVERED IN THIS CHAPTER fft ifft conj
III. DATA ANALYSIS WITH MATLAB
11 Frequency Analysis Part I: Fourier Decomposition
11.1 GOALS OF THIS CHAPTER This chapter introduces the most common method of decomposing a time series into frequency components, Fourier analysis. You will learn about the Fourier transform and the associated amplitude and phase spectra. The MATLAB® implementation of the fast Fourier transform (FFT), an efficient algorithm for calculating Fourier transformations, will be introduced and applied to the analysis of human speech sounds.
11.2 BACKGROUND Figure 11.1 shows typical recordings of two human vowel sounds. How can you characterize these different sounds? Frequency analysis provides a way to examine the relative contributions of various frequencies to an overall signal. In the case of an auditory signal, a given frequency component would be termed pitch.
11.2.1 Real Fourier Series Take some continuous function f. We can approximate such a function with a weighted series of sinusoids of various frequencies. Such a series is termed the real Fourier series: fðtÞ 5
N N X X a0 1 an cosðntÞ 1 bn sinðntÞ 2 n51 n51
ð11:1Þ
Here, the coefficients an and bn represent the relative strength of each frequency component n/2π. [a0 represents the nonoscillatory component of f(t).] So, given f(t), determining the coefficients an and bn allows for the representation of f(t) as a series sum of sinusoids.
MATLAB® for Neuroscientists. DOI: http://dx.doi.org/10.1016/B978-0-12-383836-0.00011-4
229
© 2014 Elsevier Inc. All rights reserved.
230
11. FREQUENCY ANALYSIS PART I: FOURIER DECOMPOSITION
0.2 0.15 Amplitude
0.1 0.05 0 −0.05 −0.1 −0.15
0
0.5
1
1.5
2
2.5
3
3.5
2.5
3
3.5
Time, in secs 0.15 0.1 Amplitude
0.05 0 −0.05 −0.1 −0.15 −0.2
0
0.5
1
1.5
2
Time, in secs
FIGURE 11.1
Acoustic time series representing two different human vowel sounds.
We will exploit two special properties of the sine and cosine functions to find the Fourier series coefficients an and bn. Over the interval 2 π to π, cosine and sine functions with differing frequencies have the special property of orthogonality. The integral of the product of two mutually orthogonal functions evaluates to zero. So, the integral of the product of cosine or sine functions with differing frequencies results in zero over this interval. Another interesting property of sine and cosine is that the integral of the square of a cosine or sine function over this integral is π. To find the strength, am, of a cosine component m, multiply by the corresponding cosine function and integrate: ðπ
ðπ fðtÞ cos ðmtÞdt5
2π
2π
N ð N ð X X a0 cos ðmtÞdt1 cos ðmtÞan cos ðntÞdt1 cos ðmtÞbn sin ðntÞdt 2 n51 n51 π
π
2π
2π
ð11:2Þ All terms on the right side except the cosine term where m 5 n yield zero: ðπ
ðπ fðtÞ cos ðmtÞdt 5 am 2π
cos2 ðmtÞdt 2π
III. DATA ANALYSIS WITH MATLAB
ð11:3Þ
231
11.3 EXERCISES
The right side integral evaluates to one over the integration range, yielding an expression for the Fourier series term coefficient: ðπ fðtÞ cos ðmtÞdt 5 πam
ð11:4Þ
2π
am 5
1 π
ðπ fðtÞ cos ðmtÞdt
ð11:5Þ
2π
In general, the interval of f(t) will not be 2 π to π. For an interval centered on x with length 2L, the expression becomes 1 am 5 L
x1L ð
fðtÞ cos
π L
mt dt
ð11:6Þ
x2L
A similar procedure using sine functions yields the coefficients for the sine terms of the Fourier series.
11.3 EXERCISES
EXERCISE 11.1 Write a MATLAB function to calculate coefficients for a real Fourier transform. Hint: The function will need to shift the
interval so that the interval encompasses the entire time series. In other words, x 5 0 and L 5 half the range of t.
11.3.1 Complex Fourier Transform Euler’s identity, eiωt 5 cos ωt 1 i sin ωt
ð11:7Þ
provides a straightforward way to formulate complex Fourier series representation for a given function, f(t): fðtÞ 5
N X
cn eint
ð11:8Þ
n52N
Similar to the real transform, coefficients for the complex Fourier transform can be found by
III. DATA ANALYSIS WITH MATLAB
232
11. FREQUENCY ANALYSIS PART I: FOURIER DECOMPOSITION
1 cm 5 2π
ðπ
fðtÞe2imt dt
ð11:9Þ
2π
for a given coefficient m over the interval 2 π to π. Over the interval x 2 L to x 1 L, this becomes 1 Cm 5 2L
x1L ð
fðtÞe2ðiπmt=L Þdt
ð11:10Þ
x2L
EXERCISE 11.2 Write a MATLAB function to calculate coefficients for a complex Fourier transform. This is essentially the discrete Fourier transform (DFT): Fk 5
X 1 N21 2πn fn e2i N k N n50
where k ranges from 0 to N 2 1, N is the number of points, and fn is the value of the function at point n.
ð11:11Þ
Let’s look at how this method of the Fourier transform scales with N. Given a time series with N values, this method requires a multiplication of the series and the corresponding Fourier component and subsequent sum for each coefficient. Assuming a number of coefficients equivalent to N, then you have a process that scales with N2. In other words, as N increases, the time required to compute the Fourier transform increases as N2.
11.3.2 Fast Fourier Transform With a few special tricks, a faster algorithm, the fast Fourier transform (FFT) that scales in N log N time can be formulated. One of these tricks involves taking advantage of datasets exactly 2N elements long. The increase in processing speed has made the FFT ubiquitous in signal processing. While a complete derivation of the algorithm is beyond the scope of this book, invoking the MATLAB implementation of the FFT will be discussed. MATLAB provides an FFT function fft(X), where X is a vector in time space. fft returns the frequency space representation of X. To visualize the importance of the difference in scaling, execute the following code: figure hold on N 5 1:10 * 100; plot(N, N.^ 2, 'b') plot(N, N.*log(N), 'r')
III. DATA ANALYSIS WITH MATLAB
11.3 EXERCISES
233
EXERCISE 7.3 If N represents sample size, what can you observe about the benefits of scaling as N grows? Where does the efficiency of the
FFT algorithm benefit most, for large N or small N?
11.3.3 The Inverse DFT As you might imagine, there is an inverse to the DFT: fn 5
N 21 X
2πn
Fk ei N k
ð11:12Þ
k50
MATLAB provides ifft() to perform the inverse discrete Fourier transform. EXERCISE 7.4 Generate a single sine wave. Use fft() to generate the discrete Fourier transform. Use
ifft() to retrieve the original sine wave from the DFT.
11.3.4 Amplitude Spectrum Often when you are using Fourier analysis, the amplitude spectrum is one of the first analyses performed. The amplitude spectrum graphs amplitude against frequency. In terms of the Fourier series representation, the amplitude spectrum depicts the magnitude of the coefficients at various frequencies. As such, it depicts the relative strengths of the various frequency components. The following code generates a time series composed of 10 sine waves whose frequencies and amplitudes vary systematically. L 5 1000; X 5 zeros(1,L); sampling_interval 5 0.1; t 5 (1:L) * sampling_interval; for N 5 1:10 X 5 X 1 N * sin (N*pi*t); end plot(t, X); Y 5 fft(X)/L; Now, the variable Y contains the normalized FFT of X. Note the normalization factor L. Displaying the amplitude spectrum of X requires plotting the amplitudes at various
III. DATA ANALYSIS WITH MATLAB
234
11. FREQUENCY ANALYSIS PART I: FOURIER DECOMPOSITION
frequencies. Note that fft returns only a single value, the transform coefficients. Now, how do you determine the frequency scale? The return value of the FFT assumes that frequency is evenly spaced, from 0 to a theoretical result called the Nyquist limit. Nyquist demonstrated that a discrete sampling of a continuous process can capture frequencies no higher than half the sampling frequency. Since the code above has the sampling interval, this Nyquist limit is half the inverse of the sampling interval. The following code calculates the Nyquist limit for the time series: NyLimit 5 (1 / sampling_interval)/ 2; When viewing the FFT, it is important to remember that the result is the complex transform. Thus, simply using the result of the FFT as a set of real coefficients can cause a number of problems. To display the amplitude spectrum, the absolute value of the complex coefficients will be used. The values returned by fft are the coefficients for frequencies from the negative Nyquist limit to the positive Nyquist limit. If the time series data are purely real, then the resultant transform will have even symmetry. That is, the transform will be symmetrical across the abscissa. So, in this very frequent case, only the first half of the result of fft is used. The following code employs linspace to generate frequency values and plots the amplitude spectrum. linspace generates a linearly spaced sequence of values given initial and final values. Here, the initial and final values are 0 and 1, with a value count of L/2. The resultant vector is scaled by the Nyquist limit to generate the frequency vector. F 5 linspace(0,1,L/2)*NyLimit; plot(F, abs(Y(1:L/2)));
11.3.5 Power Power at a given frequency is defined as ΦðωÞ 5 jFðωÞj2 5 FðωÞF ðωÞ
ð11:13Þ
where F* is the complex conjugate of F. To do this in MATLAB, use the function conj to return the complex conjugate of a series of complex values. Here is a plot of the power spectrum of the time series generated for the amplitude spectrum: plot(F, (Y(1:L/2).*conj(Y(1:L/2)));
11.3.6 Phase Analysis and Coherence A power spectrum alone is not a complete representation of the information in the original signal. The various Fourier components can have various phases relative to one another, as illustrated in Figure 11.2. You can plot relative phase by frequency by plotting the inverse tangent of the ratio between the imaginary component and the real component. Why is this the case? Imagine
III. DATA ANALYSIS WITH MATLAB
235
11.3 EXERCISES
FIGURE 11.2
Phase difference between two sinusoids with the same frequency.
2 sin(x) sin(x+1)
1.5 1 0.5 0 −0.5 −1 −1.5 −2
0
1
2
3
4
5
6
7
8
9
10
the complex plane, with pure real values along the abscissa (x-axis) and pure imaginary values along the ordinate (y-axis). Any complex value in your 1D Fourier transform can be represented with a coordinate pair. The magnitude of the value is simply the distance from the origin to the coordinates, or the complex modulus. The phase is the angle formed by the abscissa and the line passing through theorigin and the complex point. Thus, using basic trigonometry, the phase angle is tan21 imag real . How can you represent this in MATLAB? L 5 1000; X 5 zeros(1,L); sampling_interval 5 0.1; t 5 (1:L) * sampling_interval; for N 5 1:10 X 5 X 1 N * sin (N*pi*t); end plot(t, X); Y 5 fft(X)/L; phi 5 atan(imag(Y)./real(Y)); F 5 linspace(0,1,L/2)*NyLimit; plot(F, phi(1:L/2)); EXERCISE 11.5 Compare the phase spectrum generated in the preceding exercises with the phase
spectrum of the corresponding cosine function. Compare their power spectra.
III. DATA ANALYSIS WITH MATLAB
236
11. FREQUENCY ANALYSIS PART I: FOURIER DECOMPOSITION
TABLE 11.1 Average First and Second Formant Frequencies for Selected American English Vowels Vowel sound
First formant
Second formant
Bit
342
2322
But
623
1200
Bat
588
1952
Boot
378
997
(Data from Hillenbrand et al., 1995.)
11.4 PROJECT In this project, you will be asked to use Fourier decomposition to analyze vowel sounds produced by human speakers. On the companion website, you will find five examples of vowel sounds as produced by male American English speakers. Each sound corresponds to one of the vowel sounds in Table 11.1. The formant frequencies in Table 11.1 note the average formant frequencies as spoken by a male speaker of American English. You will use power spectra of these sounds to classify the recordings as one of these vowel sounds in the table. To complete this project, you need to understand how formants relate to frequency analysis. The human vocal tract has multiple cavities in which speech sounds resonate. As such, most sounds have multiple strong frequency components. In classifying speech sounds, the lowest strong frequency band is termed the first formant. The next highest is termed the second formant, and so on. Vowels lend themselves to a particularly simple characterization through their formants. Typically, vowel sounds have distinguishable first and second formants. Table 11.1 shows first and second formants for four vowel sounds in American English. Thus, the short “i” sound would have strong frequency representation at 342 Hz and at 2322 Hz.
MATLAB FUNCTIONS, COMMANDS, AND OPERATORS COVERED IN THIS CHAPTER fft ifft conj
III. DATA ANALYSIS WITH MATLAB