You are here : matlab → Signal Processing → pmtm
# pmtm() - Signal Processing

### Syntax

### Example

### Output / Return Value

### Limitations

### Alternatives / See Also

### Reference

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.

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

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)