You are here : matlab → Signal Processing → pmtm

pmtm() - Signal Processing

examplepxx = 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 2nw – 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.
examplepxx = 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
2nw – 1
Slepian tapers in the PSD estimate.
examplepxx = 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π).
example[pxx,f] = pmtm(___,fs) returns
a frequency vector, f, in cycles per unit time.
The sampling frequency, 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).

[pxx,w] = pmtm(x,nw,w) returns
the two-sided multitaper PSD estimates at the normalized frequencies
specified in the vector, w. w must
contain at least two elements.
[pxx,f] = pmtm(x,nw,f,fs) returns
the two-sided multitaper PSD estimates at the frequencies specified
in the vector, f. f must
contain at least two elements. The frequencies in f are
in cycles per unit time. The sampling frequency, fs,
is the number of samples per unit time. If the unit of time is seconds,
then f is in cycles/second (Hz).

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

example[___] = 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}).

example[___] = 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.
example[___] = pmtm(___,freqrange) returns
the multitaper PSD estimate over the frequency range specified by freqrange.
Valid options for freqrange are 'onesided', 'twosided',
and 'centered'.
example[___,pxxc] = pmtm(___,'ConfidenceLevel',probability) returns
the probability × 100%
confidence intervals for the PSD estimate in pxxc. 

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

Syntax

pxx = pmtm(x) examplepxx = pmtm(x,nw) examplepxx = pmtm(x,nw,nfft) example[pxx,w] = pmtm(___)[pxx,f] = pmtm(___,fs) example[pxx,w] = pmtm(x,nw,w)[pxx,f] = pmtm(x,nw,f,fs)[___] = pmtm(___,method) example[___] = pmtm(x,e,v) example[___] = pmtm(x,dpss_params)[___] = pmtm(___,'DropLastTaper',dropflag) example[___] = pmtm(___,freqrange) example[___,pxxc] = pmtm(___,'ConfidenceLevel',probability) examplepmtm(___) example

Example

Multitaper Estimate Using Default InputsOpen This Example
Obtain the multitaper PSD estimate of an input signal consisting of a discrete-time sinusoid with an angular frequency of 
 rad/sample with additive N(0,1) white noise.
Create a sine wave with an angular frequency of 
 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 ProductOpen This Example
Obtain the multitaper PSD estimate with a specified time-halfbandwidth product.
Create a sine wave with an angular frequency of 
 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 
 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 LengthOpen This Example
Obtain the multitaper PSD estimate of an input signal consisting of a discrete-time sinusoid with an angular frequency of 
 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 
 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 RateOpen This ExampleObtain 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 DFT length equal to the signal length.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 WeightsOpen This Example
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 ConcentrationsOpen This Example
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')

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

idx =

     5


Plot the multitaper PSD estimate.figure
pmtm(x,e(:,1:idx),v(1:idx),length(x),fs,'DropLastTaper',false)

DC-Centered Multitaper PSD EstimateOpen This ExampleObtain 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 BoundsOpen This Example
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 SignalOpen This Example
Generate 1024 samples of a multichannel signal consisting of three sinusoids in additive 
 white Gaussian noise. The sinusoids' frequencies are 
, 
, and 
 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)

pmtm() - Signal Processing

Syntax

Example

Output / Return Value

Limitations

Alternatives / See Also

Reference