1 % spectrum.m -
compute fft of X(4848:5360) (512 samples around 5104)
3 clear Pyy
4 W = hamming(512);
% 32 ms window at 16kHz SR
5 Y = fft(W .* X(4848:5360),512);
6 Pyy = (Y .* conj(Y)/512)';
8 S = log(Pyy(1:256));
9 f = 16000 * (0:255)/512;
10 plot(f,S(1:256))
12 xlabel('f (Hz)')
The program assumes that an input signal is in the variable, X. The output spectrum will be in the variable Pyy; line 3 clears this variable, in case it has been used previously for some other spectrum.
Line 4 defines a 512-sample Hamming window. (Fig. 4.1 shows a Hanning window, which is similar.)
In line 5, the input portion X(4848:5360) is multiplied by the Hamming window W (see fig. 4.2, from a different vowel). This has the effect of diminishing to almost zero the samples near the start and end of the selection, and keeping the samples near the middle of the window close to their original value. The operation in line 5 also applies the fft function to the result, to yield its Fourier transform, Y.
Y is a vector of complex numbers. Line 6 computes the magnitude of that vector, Pyy. The spectral magnitude is symmetrical (try plot Pyy to see this) - we are only interested in its lower half. Line 8 takes the logarithm of its lower half (samples 1:256), S, and line 9 calculates a suitable frequency scale f, i.e. in order to convert the cell number of the vector into a frequency in Hz. Line 10 then plots S against f, and line 12 labels the x axis.
If you write [f, S] to a text file, it could be processed by other software. For example, by importing it into Microsoft Excel, it is easy to plot a graph similar to figure 4.3 .