-
Notifications
You must be signed in to change notification settings - Fork 4
/
mtspec.m
153 lines (121 loc) · 3.25 KB
/
mtspec.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
function [spec,freq,sk,yk,se] = mtspec(x,dt,tbp,kspec,method,pad)
%
% Calculate the multitaper spectral estimate of time series x
%
%
% Input
%
% x(npts) - real, data vector
% dt - real, sampling rate of the time series x
% [1.0 sps - default]
% tbp - real, the time-bandwidth product (NW)
% [4.0 - default]
% kspec - integer, number of tapers to use in the estimation
% [2*tbp-1 - default]
% method - 0 - Constant; 1 - Adaptive, 2 - Quadratic
% [0 - default]
% pad - Choose if zero-padding is desired. 0 - no pad
% [0 - default]
% Output
%
% spec real vector with the adaptive estimated spectrum
% freq real vector with frequency bins
% sk eigenspectra (yk**2)
% yk eigencoefficients
% se number of degrees of freedom for each frequency bin
%
% TO DO
% Confidence intervals
% Lines, and line removal
% Other stuff
%
% Optional Input
%
% rshape Perform F-test for lines, and reshape spectrum.
% If rshape=1, then don't put the lines back
% (but keep units correct)
% fcrit The threshold probability for the F test
%
% Optional Output
%
% err jackknife 95% confidence interval
% wt array containing the ne weights for kspec eigenspectra
% normalized so that the sum of squares over
% the kspec eigenspectra is one
% fstat F statistics for single line
%------------------------------------
% Check input/Outputs
if (nargin > 6)
error('MTSPEC - Too many input arguments')
end
if (nargin < 2)
dt = 1.0;
end
if (nargin < 3)
tbp = 4.0;
end
if (nargin < 4)
kspec = floor(2*tbp)-1;
end
if (nargin < 5)
method = 0;
end
if (nargin < 6)
pad = 0;
end
if (nargout > 5)
error('Too many output arguments')
end
%------------------------------------
% Length time series and frequency
x = x(:);
npts = length(x);
nfft = npts;
if (pad~=0)
nfft = 2*npts+1;
end
[freq,nf,df] = create_fvector(x,nfft,dt);
%---------------------------------------
% Remove mean of signal
xmean = mean(x);
xvar = var(x);
x2 = detrend(x,'constant');
%----------------------------------------
% Get DPSS
[vn,lambda] = dpss(npts,tbp,kspec);
%----------------------------------------
% Get eigenspectra
[yk,sk] = eigenspec(x2,vn,kspec,nf,npts,nfft);
%----------------------------------------
% Adaptive weighted spectrum
if (method==0)
[spec,se,wt] = avgspec(nfft,nf,kspec,lambda,yk,sk);
else
[spec,se,wt] = adaptspec(nfft,nf,kspec,lambda,sk);
end
%----------------------------------------
% Jackknife spectrum
%
% !!!!!!!!! TO DO !!!!!!!!!
%
%
%----------------------------------------
% The Quadratic Inverse problem
if (method==2)
[spec,ds,dds] = qiinv(npts,nf,tbp,kspec,lambda,vn,yk,wt,spec);
end
%----------------------------------------
% Keep power relative to variance
% double power in positive frequencies
if (isreal(x))
spec(2:nf-1) = 2.d0 * spec(2:nf-1);
end
sscal = sum(spec);
sscal = xvar/(sscal*df);
spec = sscal * spec;
%-----------------------------------------
% Change to vertical vector
if (size(spec,1)<size(spec,2))
spec = spec';
end
return