在MATLAB v6.5中有spectrum的函数：
help spectrum
SPECTRUM Power spectrum estimate of one or two data sequences.
P=SPECTRUM(X,NFFT,NOVERLAP,WIND) estimates the Power Spectral Density of
signal vector X using Welch's averaged periodogram method. The signal X
is divided into overlapping sections, each of which is detrended and
windowed by the WINDOW parameter, then zero padded to length NFFT. The
magnitude squared of the length NFFT DFTs of the sections are averaged
to form Pxx.P is a two column matrix P = ; the second column
Pxxc is the 95% confidence interval. The number of rows of P is NFFT/2+1
for NFFT even, (NFFT+1)/2 for NFFT odd, or NFFT if the signal X is comp-
lex. If you specify a scalar for WINDOW, a Hanning window of that length
is used.
= SPECTRUM(X,NFFT,NOVERLAP,WINDOW,Fs) given a sampling frequency
Fs returns a vector of frequencies the same length as Pxx at which the
PSD is estimated.PLOT(F,P(:,1)) plots the power spectrum estimate
versus true frequency.
= SPECTRUM(X,NFFT,NOVERLAP,WINDOW,Fs,Pr) where Pr is a scalar
between 0 and 1, overrides the default 95% confidence interval and
returns the Pr*100% confidence interval for Pxx instead.
SPECTRUM(X) with no output arguments plots the PSD in the current
figure window, with confidence intervals.
The default values for the parameters are NFFT = 256 (or LENGTH(X),
whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2,
and Pr = .95. You can obtain a default parameter by leaving it out
or inserting an empty matrix [], e.g. SPECTRUM(X,[],128).
P = SPECTRUM(X,Y) performs spectral analysis of the two sequences
X and Y using the Welch method. SPECTRUM returns the 8 column array
P =
where
Pxx= X-vector power spectral density
Pyy= Y-vector power spectral density
Pxy= Cross spectral density
Txy= Complex transfer function from X to Y = Pxy./Pxx
Cxy= Coherence function between X and Y = (abs(Pxy).^2)./(Pxx.*Pyy)
Pxxc,Pyyc,Pxyc = Confidence range.
All input and output options are otherwise exactly the same as for the
single input case.
SPECTRUM(X,Y) with no output arguments will plot Pxx, Pyy, abs(Txy),
angle(Txy) and Cxy in sequence, pausing between plots.
SPECTRUM(X,...,DFLAG), where DFLAG can be 'linear', 'mean' or 'none',
specifies a detrending mode for the prewindowed sections of X (and Y).
DFLAG can take the place of any parameter in the parameter list
(besides X) as long as it is last, e.g. SPECTRUM(X,'none');
See also PSD, CSD, TFE, COHERE, SPECGRAM, SPECPLOT, DETREND, PMTM,
PMUSIC.
ETFE, SPA, and ARX in the Identification Toolbox.
很明显，文中主要是用psdfun函数进行功率谱分析，并提取信号的特征:
function y=psdfun(x)%函数用来求信号的特征频率向量
=spectrum(x,512,128,hamming(216),22050);
m=0;
for k=2:length(p)-1;
if p(k-1)<p(k)&p(k)>p(k+1)
m=m+1;
y(m)=f(k);
end
end
从上贴出的spectrum函数可知，该函数是用Welch方法求平均功率谱，函数中x是信号，变换长度NFFT=512，非重叠长度noverlap=128，窗函数hamming(216)及采样频率Fs=22050。然后用p计算局部峰值，存贮在数组y中。
正因为是求平均功率谱，使频谱比较光滑，才有5~6页中的谱图。