Exercise 4 (1 Point):

Multitaper single trial confidence intervals.

Generate a single (4-second long) realization of a f = 40 Hz sinusoid plus Gaussian white noise (zero mean, unit variance) X(t)= cos(2πft)+ε(t). Use a sampling rate of 1000 Hz.

R: See code below

Compute the multitaper power spectral function estimate using a bandwidth of 4 Hz. Plot also the corresponding 95% confidence intervals based on the $\chi^2$ approximation.

R: See code below

Based on Lecture 4, how does the performance of the $\chi^2$ approximation depends on the number of degrees of freedom?

R: For $\nu > 20$ degrees of freedom the $\chi^2$ approximation agrees with bootstrap approximation



In [1]:

    
using PlotlyJS
using DSP
using StatsFuns  #for inverse of chi-cdf









    






 
 Plotly javascript loaded.
 To load again call 
init_notebook(true)



In [2]:

    
# Define signal
f = 40 
Fs = 1000
dt = 1/Fs
T = 4
df = 1/T
N = T*Fs
fNyquist= Fs/2
t = collect(0:N-1)*dt
s = cos.(2*π*f*t) .+ randn(N);



In [3]:

    
# Multitaper   
# Note: mt_pgram doesn't take care of removing sample mean and temporal mean as in MATLAB. 
# Do we need to do it manually?
R = 4 #bandwidth
ntapers = R*T - 1
Smt = power(mt_pgram(s; fs=Fs, nw=R, ntapers = ntapers))
Smt_db = 10*log10.(Smt/maximum(Smt));



In [4]:

    
F = collect(0:df:fNyquist) #frequency grid
#confidence intervals
ntrials = 1
dof = 2*ntapers*ntrials
alpha=0.05
q1=chisqinvcdf(dof, alpha/2);
q2=chisqinvcdf(dof, 1-alpha/2);
uci=dof.*Smt_db./q2;
lci=dof.*Smt_db./q1;

lci_trace = PlotlyJS.scatter(; x = F, y= lci, marker_color="blue", line_width=0.2, name="Lower CI")
trace = PlotlyJS.scatter(x = F, y=Smt_db)
uci_trace = PlotlyJS.scatter(; x = F, y= uci, marker_color="blue", line_width=0.2, name="Upper CI")


layout = Layout(;title="Multitaper", showlegend=false,
                         xaxis=attr(title="Frequency (Hz)", zeroline=false, range=[20, 60]),
                         yaxis=attr(title="Power (dB)", zeroline=false))

data = [lci_trace, trace, uci_trace]
p = plot(data, layout)









    Out[4]: