Documentation

pxx = pmtm(x) returns Thomson’s multitaper power spectral density (PSD) estimate, pxx, of the input signal, x. When x is a vector, it is treated as a single channel. When x is a matrix, the PSD is computed independently for each column and stored in the corresponding column of pxx. The tapers are the discrete prolate spheroidal (DPSS), or Slepian, sequences. The time-halfbandwidth, nw, product is 4. By default, pmtm uses the first 2 × nw – 1 DPSS sequences. If x is real-valued, pxx is a one-sided PSD estimate. If x is complex-valued, pxx is a two-sided PSD estimate. The number of points, nfft, in the discrete Fourier transform (DFT) is the maximum of 256 or the next power of two greater than the signal length.

pxx = pmtm(x,nw) use the time-halfbandwidth product, nw, to obtain the multitaper PSD estimate. The time-halfbandwidth product controls the frequency resolution of the multitaper estimate. pmtm uses 2 × nw – 1 Slepian tapers in the PSD estimate.

pxx = pmtm(x,nw,nfft) uses nfft points in the DFT. If nfft is greater than the signal length, x is zero-padded to length nfft. If nfft is less than the signal length, the signal is wrapped modulo nfft.

[pxx,w] = pmtm(___) returns the normalized frequency vector, w. If pxx is a one-sided PSD estimate, w spans the interval [0,π] if nfft is even and [0,π) if nfft is odd. If pxx is a two-sided PSD estimate, w spans the interval [0,2π).

[pxx,f] = pmtm(___,fs) returns a frequency vector, f, in cycles per unit time. The sample rate, fs, is the number of samples per unit time. If the unit of time is seconds, then f is in cycles/sec (Hz). For real–valued signals, f spans the interval [0,fs/2] when nfft is even and [0,fs/2) when nfft is odd. For complex-valued signals, f spans the interval [0,fs). fs must be the fourth input to pmtm. To input a sample rate and still use the default values of the preceding optional arguments, specify these arguments as empty, [].

[pxx,w] = pmtm(x,nw,w) returns the two-sided multitaper PSD estimates at the normalized frequencies specified in w. The vector w must contain at least two elements, because otherwise the function interprets it as nfft.

[pxx,f] = pmtm(x,nw,f,fs) returns the two-sided multitaper PSD estimates at the frequencies specified in the vector, f. The vector f must contain at least two elements, because otherwise the function interprets it as nfft. The frequencies in f are in cycles per unit time. The sample rate, fs, is the number of samples per unit time. If the unit of time is seconds, then f is in cycles/second (Hz).

[___] = pmtm(___,method) combines the individual tapered PSD estimates using the method, method. method can be one of: 'adapt' (default), 'eigen', or 'unity'.

[___] = pmtm(x,e,v) uses the tapers in the N-by-K matrix e with concentrations v in the frequency band [-w,w]. N is the length of the input signal, x. Use dpss to obtain the Slepian tapers and corresponding concentrations.

[___] = pmtm(x,dpss_params) uses the cell array, dpss_params, to pass input arguments to dpss except the number of elements in the sequences. The number of elements in the sequences is the first input argument to dpss and is not included in dpss_params. An example of this usage is pxx = pmtm(randn(1000,1),{2.5,3}).

[___] = pmtm(___,'DropLastTaper',dropflag) specifies whether pmtm drops the last taper in the computation of the multitaper PSD estimate. dropflag is a logical. The default value of dropflag is true and the last taper is not used in the PSD estimate.

[___] = pmtm(___,freqrange) returns the multitaper PSD estimate over the frequency range specified by freqrange. Valid options for freqrange are 'onesided', 'twosided', and 'centered'.

[___,pxxc] = pmtm(___,'ConfidenceLevel',probability) returns the probability × 100% confidence intervals for the PSD estimate in pxxc.

pmtm(___) with no output arguments plots the multitaper PSD estimate in the current figure window.

Examples

Multitaper Estimate Using Default Inputs

Obtain the multitaper PSD estimate of an input signal consisting of a discrete-time sinusoid with an angular frequency of $π / 4$ rad/sample with additive N(0,1) white noise.

Create a sine wave with an angular frequency of $π / 4$ rad/sample with additive N(0,1) white noise. The signal is 320 samples in length. Obtain the multitaper PSD estimate using the default time-halfbandwidth product of 4 and DFT length. The default number of DFT points is 512. Because the signal is real-valued, the PSD estimate is one-sided and there are 512/2+1 points in the PSD estimate.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
pxx = pmtm(x);

Plot the multitaper PSD estimate.

pmtm(x)

Specify Time-Halfbandwidth Product

Obtain the multitaper PSD estimate with a specified time-halfbandwidth product.

Create a sine wave with an angular frequency of $π / 4$ rad/sample with additive N(0,1) white noise. The signal is 320 samples in length. Obtain the multitaper PSD estimate with a time-halfbandwidth product of 2.5. The resolution bandwidth is $[- 2.5 π / 320, 2.5 π / 320]$ rad/sample. The default number of DFT points is 512. Because the signal is real-valued, the PSD estimate is one-sided and there are 512/2+1 points in the PSD estimate.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
pmtm(x,2.5)

DFT Length Equal to Signal Length

Obtain the multitaper PSD estimate of an input signal consisting of a discrete-time sinusoid with an angular frequency of $π / 4$ rad/sample with additive N(0,1) white noise. Use a DFT length equal to the signal length.

Create a sine wave with an angular frequency of $π / 4$ rad/sample with additive N(0,1) white noise. The signal is 320 samples in length. Obtain the multitaper PSD estimate with a time-halfbandwidth product of 3 and a DFT length equal to the signal length. Because the signal is real-valued, the one-sided PSD estimate is returned by default with a length equal to 320/2+1.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
pmtm(x,3,length(x))

Multitaper Estimate with Sample Rate

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+randn(size(t));
[pxx,f] = pmtm(x,3,length(x),fs);

Plot the multitaper PSD estimate.

pmtm(x,3,length(x),fs)

Average Single-Taper Estimates with Unity Weights

Obtain a multitaper PSD estimate where the individual tapered direct spectral estimates are given equal weight in the average.

Obtain the multitaper PSD estimate of a signal sampled at 1 kHz. The signal is a 100 Hz sine wave in additive N(0,1) white noise. The signal duration is 2 s. Use a time-halfbandwidth product of 3 and a DFT length equal to the signal length. Use the 'unity' option to give equal weight in the average to each of the individual tapered direct spectral estimates.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+randn(size(t));
[pxx,f] = pmtm(x,3,length(x),fs,'unity');

Plot the multitaper PSD estimate.

pmtm(x,3,length(x),fs,'unity')

DPSS Sequences and Their Frequency-Domain Concentrations

This example examines the frequency-domain concentrations of the DPSS sequences. The example produces a multitaper PSD estimate of an input signal by precomputing the Slepian sequences and selecting only those with more than 99% of their energy concentrated in the resolution bandwidth.

The signal is a 100 Hz sine wave in additive N(0,1) white noise. The signal duration is 2 s.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+randn(size(t));

Set the time-halfbandwidth product to 3.5. For the signal length of 2000 samples and a sampling interval of 0.001 seconds, this results in a resolution bandwidth of [-1.75,1.75] Hz. Calculate the first 10 Slepian sequences and examine their frequency concentrations in the specified resolution bandwidth.

[e,v] = dpss(length(x),3.5,10);

stem(1:length(v),v,'filled')
ylim([0 1.2])
title('Proportion of Energy in [-w,w] of k-th Slepian Sequence')

Determine the number of Slepian sequences with energy concentrations greater than 99%. Using the selected DPSS sequences, obtain the multitaper PSD estimate. Set 'DropLastTaper' to false to use all the selected tapers.

hold on
plot(1:length(v),0.99*ones(length(v),1))

idx = find(v>0.99,1,'last')

idx = 5

[pxx,f] = pmtm(x,e(:,1:idx),v(1:idx),length(x),fs,'DropLastTaper',false);

Plot the multitaper PSD estimate.

figure
pmtm(x,e(:,1:idx),v(1:idx),length(x),fs,'DropLastTaper',false)

DC-Centered Multitaper PSD Estimate

Obtain the multitaper PSD estimate of a 100 Hz sine wave in additive N(0,1) noise. The data are sampled at 1 kHz. Use the 'centered' option to obtain the DC-centered PSD.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+randn(size(t));
[pxx,f] = pmtm(x,3.5,length(x),fs,'centered');

Plot the DC-centered PSD estimate.

pmtm(x,3.5,length(x),fs,'centered')

Upper and Lower 95%-Confidence Bounds

The following example illustrates the use of confidence bounds with the multitaper PSD estimate. While not a necessary condition for statistical significance, frequencies in the multitaper PSD estimate where the lower confidence bound exceeds the upper confidence bound for surrounding PSD estimates clearly indicate significant oscillations in the time series.

Create a signal consisting of the superposition of 100-Hz and 150-Hz sine waves in additive white N(0,1) noise. The amplitude of the two sine waves is 1. The sampling frequency is 1 kHz. The signal is 2 s in duration.

fs = 1000;
t = 0:1/fs:2-1/fs;
x = cos(2*pi*100*t)+cos(2*pi*150*t)+randn(size(t));

Obtain the multitaper PSD estimate with 95%-confidence bounds. Plot the PSD estimate along with the confidence interval and zoom in on the frequency region of interest near 100 and 150 Hz.

[pxx,f,pxxc] = pmtm(x,3.5,length(x),fs,'ConfidenceLevel',0.95);

plot(f,10*log10(pxx))
hold on
plot(f,10*log10(pxxc),'r-.')
xlim([85 175])
xlabel('Hz')
ylabel('dB')
title('Multitaper PSD Estimate with 95%-Confidence Bounds')

The lower confidence bound in the immediate vicinity of 100 and 150 Hz is significantly above the upper confidence bound outside the vicinity of 100 and 150 Hz.

Multitaper PSD Estimate of a Multichannel Signal

Generate 1024 samples of a multichannel signal consisting of three sinusoids in additive $N (0, 1)$ white Gaussian noise. The sinusoids' frequencies are $π / 2$ , $π / 3$ , and $π / 4$ rad/sample. Estimate the PSD of the signal using Thomson's multitaper method and plot it.

N = 1024;
n = 0:N-1;

w = pi./[2;3;4];
x = cos(w*n)' + randn(length(n),3);

pmtm(x)

Input Arguments

`x` — Input signal
vector | matrix

Input signal, specified as a row or column vector, or as a matrix. If x is a matrix, then its columns are treated as independent channels.

Example: cos(pi/4*(0:159))+randn(1,160) is a single-channel row-vector signal.

Example: cos(pi./[4;2]*(0:159))'+randn(160,2) is a two-channel signal.

Data Types: single | double
Complex Number Support: Yes

`nw` — Time-halfbandwidth product
4 (default) | positive scalar

Time-halfbandwidth product, specified as a positive scalar. In multitaper spectral estimation, the user specifies the resolution bandwidth of the multitaper estimate [–W,W] where W = k/NΔt for some small k > 1. Equivalently, W is some small multiple of the frequency resolution of the DFT. The time-halfbandwidth product is the product of the resolution halfbandwidth and the number of samples in the input signal, N. The number of Slepian tapers whose Fourier transforms are well-concentrated in [–W,W] (eigenvalues close to unity) is 2NW – 1.

`nfft` — Number of DFT points
`max(256,2^nextpow2(length(x)))` (default) | integer | `[]`

Number of DFT points, specified as a positive integer. For a real-valued input signal, x, the PSD estimate, pxx has length (nfft/2 + 1) if nfft is even, and (nfft + 1)/2 if nfft is odd. For a complex-valued input signal,x, the PSD estimate always has length nfft. If nfft is specified as empty, the default nfft is used.

Data Types: single | double

`fs` — Sample rate
positive scalar

Sample rate, specified as a positive scalar. The sample rate is the number of samples per unit time. If the unit of time is seconds, then the sample rate has units of Hz.

`w` — Normalized frequencies
vector

Normalized frequencies, specified as a row or column vector with at least two elements. Normalized frequencies are in rad/sample.

Example: w = [pi/4 pi/2]

Data Types: double

`f` — Frequencies
vector

Frequencies, specified as a row or column vector with at least two elements. The frequencies are in cycles per unit time. The unit time is specified by the sample rate, fs. If fs has units of samples/second, then f has units of Hz.

Example: fs = 1000; f = [100 200]

Data Types: double

`method` — Weights on individual tapered PSD estimates
`'adapt'` (default) | `'eigen'` | `'unity'`

Weights on individual tapered PSD estimates, specified as one of 'adapt', 'eigen', or 'unity'. The default is Thomson’s adaptive frequency-dependent weights, 'adapt'. The calculation of these weights is detailed on pp. 368–370 in [1]. The 'eigen' method weights each tapered PSD estimate by the eigenvalue (frequency concentration) of the corresponding Slepian taper. The 'unity' method weights each tapered PSD estimate equally.

`e` — DPSS (Slepian) sequences
matrix

DPSS (Slepian) sequences, specified as a N-by-K matrix where N is the length of the input signal, x. The matrix e is the output of dpss.

`v` — Eigenvalues for DPSS (Slepian) sequences
vector

Eigenvalues for DPSS (Slepian) sequences, specified as a column vector. The eigenvalues for the DPSS sequences indicate the proportion of the sequence energy concentrated in the resolution bandwidth, [-W, W]. The eigenvalues range lie in the interval (0,1) and generally the first 2NW-1 eigenvalues are close to 1 and then decrease toward 0.

`dpss_params` — Input arguments for `dpss`
cell array

Input arguments for dpss, specified as a cell array. The first input argument to dpss is the length of the DPSS sequences and is omitted from dpss_params. The length of the DPSS sequences is obtained from the length of the input signal, x.

Example: {3.5,5}

`dropflag` — Flag indicating whether to drop or keep the last DPSS sequence
`true` (default) | `false`

Flag indicating whether to drop or keep the last DPSS sequence, specified as a logical. The default is true and pmtm drops the last taper. In a multitaper estimate, the first 2NW – 1 DPSS sequences have eigenvalues close to unity. If you use less than 2NW – 1 sequences, it is likely that all the tapers have eigenvalues close to 1 and you can specify dropflag as false to keep the last taper.

`freqrange` — Frequency range for PSD estimate
`'onesided'` | `'twosided'` | `'centered'`

Frequency range for the PSD estimate, specified as a one of 'onesided', 'twosided', or 'centered'. The default is 'onesided' for real-valued signals and 'twosided' for complex-valued signals. The frequency ranges corresponding to each option are

'onesided' — returns the one-sided PSD estimate of a real-valued input signal, x. If nfft is even, pxx has length nfft/2 + 1 and is computed over the interval [0,π] rad/sample. If nfft is odd, the length of pxx is (nfft + 1)/2 and the interval is [0,π) rad/sample. When fs is optionally specified, the corresponding intervals are [0,fs/2] cycles/unit time and [0,fs/2) cycles/unit time for even and odd length nfft respectively.
'twosided' — returns the two-sided PSD estimate for either the real-valued or complex-valued input, x. In this case, pxx has length nfft and is computed over the interval [0,2π) rad/sample. When fs is optionally specified, the interval is [0,fs) cycles/unit time.
'centered' — returns the centered two-sided PSD estimate for either the real-valued or complex-valued input, x. In this case, pxx has length nfft and is computed over the interval (–π,π] rad/sample for even length nfft and (–π,π) rad/sample for odd length nfft. When fs is optionally specified, the corresponding intervals are (–fs/2, fs/2] cycles/unit time and (–fs/2, fs/2) cycles/unit time for even and odd length nfft respectively.

`probability` — Confidence interval for PSD estimate
0.95 (default) | scalar in the range (0,1)

Coverage probability for the true PSD, specified as a scalar in the range (0,1). The output, pxxc, contains the lower and upper bounds of the probability × 100% interval estimate for the true PSD.

Output Arguments

`pxx` — PSD estimate
vector | matrix

PSD estimate, returned as a real-valued, nonnegative column vector or matrix. Each column of pxx is the PSD estimate of the corresponding column of x. The units of the PSD estimate are in squared magnitude units of the time series data per unit frequency. For example, if the input data is in volts, the PSD estimate is in units of squared volts per unit frequency. For a time series in volts, if you assume a resistance of 1 Ω and specify the sample rate in hertz, the PSD estimate is in watts per hertz.

Data Types: single | double

`w` — Normalized frequencies
vector

Normalized frequencies, returned as a real-valued column vector. If pxx is a one-sided PSD estimate, w spans the interval [0,π] if nfft is even and [0,π) if nfft is odd. If pxx is a two-sided PSD estimate, w spans the interval [0,2π). For a DC-centered PSD estimate, w spans the interval (–π,π] for even nfft and (–π,π) for odd nfft.

Data Types: double

`f` — Cyclical frequencies
vector

Cyclical frequencies, returned as a real-valued column vector. For a one-sided PSD estimate, f spans the interval [0,fs/2] when nfft is even and [0,fs/2) when nfft is odd. For a two-sided PSD estimate, f spans the interval [0,fs). For a DC-centered PSD estimate, f spans the interval (–fs/2, fs/2] cycles/unit time for even length nfft and (–fs/2, fs/2) cycles/unit time for odd length nfft.

Data Types: double | single

`pxxc` — Confidence bounds
matrix

Confidence bounds, returned as a matrix with real-valued elements. The row size of the matrix is equal to the length of the PSD estimate, pxx. pxxc has twice as many columns as pxx. Odd-numbered columns contain the lower bounds of the confidence intervals, and even-numbered columns contain the upper bounds. Thus, pxxc(m,2*n-1) is the lower confidence bound and pxxc(m,2*n) is the upper confidence bound corresponding to the estimate pxx(m,n). The coverage probability of the confidence intervals is determined by the value of the probability input.

Data Types: single | double

More About