
% clear
close all

% Matlab program to compare classical "periodogram" with the multitaper
% approach
% Examples follows Matlab doc (periodogram, pmtm)

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); 
y = x + 2*randn(size(t));     % Sinusoids plus noise

figure
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')


% now estimate the amplitude spectra

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);

% Plot single-sided amplitude spectrum.
figure
plot(f,2*abs(Y(1:NFFT/2))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

% .. .and the power spectra (simply the square
figure
plot(f,4*abs(Y(1:NFFT/2)).^2) 
title('Single-Sided Power Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

% .. .and the power spectra (simply the square
figure
semilogy(f,4*abs(Y(1:NFFT/2)).^2) 
title('Single-Sided Power Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')


% Now the Multitaper approach
% please refer to the Matlab doc for info on the parameters of pmtm -> help
% pmtm
figure
[Pxx,Pxxc,f] = pmtm(x,3.5,512,Fs,0.99);
hpsd = dspdata.psd([Pxx Pxxc],'Fs',Fs);
plot(hpsd)

