www.gusucode.com > signal 案例源码程序 matlab代码 > signal/MultitaperMethodExample.m

    %% Multitaper Method
% The periodogram can be interpreted as filtering a length $L$ signal,
% $x_L(n)$, through a filter bank (a set of filters in parallel) of $L$ FIR
% bandpass filters. The 3 dB bandwidth of each of these bandpass filters
% can be shown to be approximately equal to $f_s/L$. The magnitude response
% of each one of these bandpass filters resembles that of a rectangular
% window. The periodogram can thus be viewed as a computation of the power
% of each filtered signal (i.e., the output of each bandpass filter) that
% uses just one sample of each filtered signal and assumes that the PSD of
% $x_L(n)$ is constant over the bandwidth of each bandpass filter.
%
% As the length of the signal increases, the bandwidth of each bandpass
% filter decreases, making it a more selective filter, and improving the
% approximation of constant PSD over the bandwidth of the filter. This
% provides another interpretation of why the PSD estimate of the
% periodogram improves as the length of the signal increases. However,
% there are two factors apparent from this standpoint that compromise the
% accuracy of the periodogram estimate. First, the rectangular window
% yields a poor bandpass filter. Second, the computation of the power at
% the output of each bandpass filter relies on a single sample of the
% output signal, producing a very crude approximation.
%
% Welch's method can be given a similar interpretation in terms of a filter
% bank. In Welch's implementation, several samples are used to compute the
% output power, resulting in reduced variance of the estimate. On the other
% hand, the bandwidth of each bandpass filter is larger than that
% corresponding to the periodogram method, which results in a loss of
% resolution. The filter bank model thus provides a new interpretation of
% the compromise between variance and resolution.
%
% Thompson's _multitaper method_ (MTM) builds on these results to provide
% an improved PSD estimate. Instead of using bandpass filters that are
% essentially rectangular windows (as in the periodogram method), the MTM
% method uses a bank of optimal bandpass filters to compute the estimate.
% These optimal FIR filters are derived from a set of sequences known as
% discrete prolate spheroidal sequences (DPSSs, also known as _Slepian
% sequences_).
%
% In addition, the MTM method provides a time-bandwidth parameter with
% which to balance the variance and resolution. This parameter is given by
% the time-bandwidth product, $NW$ and it is directly related to the number
% of tapers used to compute the spectrum. There are always $2NW-1$ tapers
% used to form the estimate. This means that, as $NW$ increases, there are
% more estimates of the power spectrum, and the variance of the estimate
% decreases. However, the bandwidth of each taper is also proportional to
% $NW$, so as $NW$ increases, each estimate exhibits more spectral leakage
% (i.e., wider peaks) and the overall spectral estimate is more biased. For
% each data set, there is usually a value for $NW$ that allows an optimal
% trade-off between bias and variance.
%
% The Signal Processing Toolbox(TM) function that implements the MTM method
% is |pmtm|. Use |pmtm| to compute the PSD of a signal.

% Copyright 2015 The MathWorks, Inc.


%%

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes
f = [150;140];            % Sinusoid frequencies

xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
pmtm(xn,4,[],fs)

%%
% By lowering the time-bandwidth product, you can increase the resolution
% at the expense of larger variance.

pmtm(xn,1.5,[],fs)

%%
% This method is more computationally expensive than Welch's method due to
% the cost of computing the discrete prolate spheroidal sequences. For long
% data series (10,000 points or more), it is useful to compute the DPSSs
% once and save them in a MAT-file. |dpsssave|, |dpssload|, |dpssdir|, and
% |dpssclear| are provided to keep a database of saved DPSSs in the
% MAT-file |dpss.mat|.