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

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);
lv = length(v);

stem(1:lv,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 lv],0.99*[1 1])
hold off

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

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

%
% Plot the multitaper PSD estimate.

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

% Copyright 2015 The MathWorks, Inc.