% eda09_21
% MultiTaper analysis

% standard time/frequency set up
N=64;
Nf=N/2+1;
Dt=1;
t=Dt*[0:N-1]';
tmax = Dt*(N-1);
fmax=1/(2*Dt);
Df=fmax/(N/2);
Dw=2*pi*Df;
f=Df*[0:N/2,-N/2+1:-1]'; 

% choose bandwidth to be +/-3 Dw
f0 = 2*Df;
w0 = 2*pi*f0;

% matrix, M, of sinc functions
M = zeros(N,N);
for m = [1:N]
for n = [1:N]
    M(m,n) = 2*w0*sinc( w0*(m-n)*Dt / pi);
end
end

% solve algebric eigenvalue problem
% eig sorts in ascending order
[V,D] = eig(M);
lambda=diag(D);
disp(sprintf('6 largest eigenvalues:'));
for m = [N-5:N]
    disp(sprintf('    %d: %f', m, lambda(m) ));
end

% four largest eigenvalues
w1 = V(:,N);
w2 = V(:,N-1);
w3 = V(:,N-2);
w4 = V(:,N-3);

% power spectral density of window functions
s1 = ((Dt^2)/tmax) * abs(fft(w1)).^2;
s2 = ((Dt^2)/tmax) * abs(fft(w2)).^2;
s3 = ((Dt^2)/tmax) * abs(fft(w3)).^2;
s4 = ((Dt^2)/tmax) * abs(fft(w4)).^2;

% plot window functions
figure(1);
clf;

winmin = -0.25;
winmax = -winmin;
subplot(4,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, winmin, winmax ]);
plot( t, w1, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w1(t)');

subplot(4,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, winmin, winmax ]);
plot( t, w2, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w2(t)');

subplot(4,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, winmin, winmax ]);
plot( t, w3, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w3(t)');

subplot(4,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, winmin, winmax ]);
plot( t, w4, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w4(t)');

% plot amplitude spectral density of window functions
figure(2);
clf;

subplot(4,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(s1))] );
plot( f(1:Nf), sqrt(s1(1:Nf)), 'k-', 'LineWidth', 2);
plot( [f0, f0], [0, max(sqrt(s1))], 'k:', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|w1(f)|');

subplot(4,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(s2))] );
plot( f(1:Nf), sqrt(s2(1:Nf)), 'k-', 'LineWidth', 2);
plot( [f0, f0], [0, max(sqrt(s2))], 'k:', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|w2(f)|');

subplot(4,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(s3))] );
plot( f(1:Nf), sqrt(s3(1:Nf)), 'k-', 'LineWidth', 2);
plot( [f0, f0], [0, max(sqrt(s3))], 'k:', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|w3(f)|');

subplot(4,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(s4))] );
plot( f(1:Nf), sqrt(s4(1:Nf)), 'k-', 'LineWidth', 2);
plot( [f0, f0], [0, max(sqrt(s4))], 'k:', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|w1(f)|');

% new, longer time interval, with same sampling interval
N2=512;
Nf2=N2/2+1;
t2=Dt*[0:N2-1]';
tmax2 = Dt*(N2-1);
fmax=1/(2*Dt);
Df2=fmax/(N2/2);
Dw2=2*pi*Df2;
f2=Df2*[0:N2/2,-N2/2+1:-1]';

% create time series and window it
ww= 2*pi*fmax/3.3;
d=cos(ww*t2);
dw0 = [[1:N], 0*[1:N2-N] ]'.*d;
dw1 = [w1', 0*[1:N2-N] ]'.*d;
dw2 = [w2', 0*[1:N2-N] ]'.*d;
dw3 = [w3', 0*[1:N2-N] ]'.*d;
dw4 = [w4', 0*[1:N2-N] ]'.*d;

% plot windowed time series
figure(3);
clf;

subplot(5,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(d), max(d) ]);
plot( t2, d, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('d(t)');

subplot(5,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(d), max(d) ]);
plot( t2, dw0, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('d(t)');

subplot(5,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(dw1), max(dw1) ]);
plot( t2, dw1, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w1(t)d(t)');

subplot(5,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(dw2), max(dw2) ]);
plot( t2, dw2, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w2(t)d(t)');

subplot(5,1,5);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax2, min(dw3), max(dw3) ]);
plot( t2, dw3, 'k-', 'LineWidth', 2 );
xlabel('time t, s');
ylabel('w3(t)d(t)');

% power spectral density of time series
ds = ((Dt^2)/tmax) * abs(fft(d)).^2;

% power spectral density of windowed time series
dw0s = ((Dt^2)/tmax) * abs(fft(dw0)).^2;
dw1s = ((Dt^2)/tmax) * abs(fft(dw1)).^2;
dw2s = ((Dt^2)/tmax) * abs(fft(dw2)).^2;
dw3s = ((Dt^2)/tmax) * abs(fft(dw3)).^2;
dw4s = ((Dt^2)/tmax) * abs(fft(dw4)).^2;

% average power spectral densities
dws = (1/3)*( dw1s + dw2s + dw3s  );

% plot power spectral density of timeseries
figure(4);
clf;

subplot(6,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(ds))] );
plot( f2(1:Nf2), sqrt(ds(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d(f)|');

subplot(6,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dw0s))] );
plot( f2(1:Nf2), sqrt(dw0s(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d0(f)|');

subplot(6,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dw1s))] );
plot( f2(1:Nf2), sqrt(dw1s(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d1(f)|');

subplot(6,1,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dw2s))] );
plot( f2(1:Nf2), sqrt(dw2s(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d2(f)|');

subplot(6,1,5);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dw3s))] );
plot( f2(1:Nf2), sqrt(dw3s(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d3(f)|');

subplot(6,1,6);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, max(sqrt(dws))] );
plot( f2(1:Nf2), sqrt(dws(1:Nf2)), 'k-', 'LineWidth', 2);
xlabel('frequency f, Hz');
ylabel('|d(f)|avg');
