-
Notifications
You must be signed in to change notification settings - Fork 8
/
spm_eeg_specest_morlet.m
126 lines (106 loc) · 4.05 KB
/
spm_eeg_specest_morlet.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
function res = spm_eeg_specest_morlet(S, data, time)
% Plugin for spm_eeg_tf implementing Morlet wavelet transform
% FORMAT res = spm_eeg_specest_morlet(S, data, time)
%
% S - input structure
% fields of S:
% S.subsample - factor by which to subsample the time axis (default - 1)
% either
% S.ncycles - Morlet wavelet factor (default - 7)
% or
% S.timeres - Fixed time window length in ms
%
% S.frequencies - vector of frequencies (default - 0-48) at optimal frequency bins
%
% Output:
% res -
% If no input is provided the plugin returns a cfg branch for itself
%
% If input is provided:
% res.fourier - the complex output of wavelet transform
% res.time - time axis
% res.freq - frequency axis
%______________________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Vladimir Litvak based on the code from Stefan Kiebel and Will Penny
% $Id: spm_eeg_specest_morlet.m 3749 2010-03-04 13:17:17Z vladimir $
%-This part if for creating a config branch that plugs into spm_cfg_eeg_tf
% Any parameters can be specified and they are then passed to the plugin
% when it's called.
%--------------------------------------------------------------------------
if nargin == 0
subsample = cfg_entry;
subsample.tag = 'subsample';
subsample.name = 'Subsample';
subsample.strtype = 'n';
subsample.num = [1 1];
subsample.val = {1};
subsample.help = {'Set to N to subsample the time axis to every Nth sample (to reduce the dataset size).'};
ncycles = cfg_entry;
ncycles.tag = 'ncycles';
ncycles.name = 'Number of wavelet cycles';
ncycles.strtype = 'n';
ncycles.num = [1 1];
ncycles.val = {7};
ncycles.help = {'Number of wavelet cycles (a.k.a. Morlet wavelet factor)',...
'This parameter controls the time-frequency trade-off',...
'Increasing it increases the frequency resolution at the expense of time resolution.'};
timeres = cfg_entry;
timeres.tag = 'timeres';
timeres.name = 'Fixed time window length';
timeres.strtype = 'r';
timeres.num = [1 1];
timeres.val = {0};
timeres.help = {'Fixed time window for all frequencies.',...
'Specify time window length in ms.',...
'Default valued of 0 specifies variable time window length'};
morlet = cfg_branch;
morlet.tag = 'morlet';
morlet.name = 'Morlet wavelet transform';
morlet.val = {ncycles, timeres, subsample};
res = morlet;
return
elseif nargin < 3
error('Three input arguments are required');
end
%-Defaults
%--------------------------------------------------------------------------
if ~isfield(S, 'subsample')
S.subsample = 1;
end
if ~isfield(S, 'ncycles')
S.ncycles = 7;
end
dt = time(end) - time(1);
if ~isfield(S, 'frequencies') || isempty(S.frequencies)
S.frequencies = (1/dt):max(1/dt, floor(dt)/dt):48;
end
%-Generate wavelets
%--------------------------------------------------------------------------
if ~isfield(S, 'timeres') || (S.timeres == 0)
M = spm_eeg_morlet(S.ncycles, 1000*diff(time(1:2)), S.frequencies);
else
M = spm_eeg_morlet(S.ncycles, 1000*diff(time(1:2)), S.frequencies, 1000./S.timeres);
end
%-Data dimensions
%--------------------------------------------------------------------------
Nchannels = size(data, 1);
Nsamples = size(data, 2);
Nfrequencies = numel(M);
%-Initialize output struct
%--------------------------------------------------------------------------
res = [];
res.freq = S.frequencies;
res.time = time(1:S.subsample:end);
res.fourier = zeros(Nchannels, Nfrequencies, length(res.time));
%-Compute wavelet transform
%--------------------------------------------------------------------------
for j = 1:Nchannels
for i = 1:Nfrequencies
tmp = conv(data(j, :), M{i});
% time shift to remove delay
tmp = tmp([1:Nsamples] + (length(M{i})-1)/2);
tmp = tmp(1:S.subsample:end);
res.fourier(j, i, :) = tmp;
end
end