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|.