-
Notifications
You must be signed in to change notification settings - Fork 0
/
rd_simSampledSinewave.m
157 lines (134 loc) · 3.72 KB
/
rd_simSampledSinewave.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
154
155
156
% rd_simSampledSinewave.m
%% setup
sinewave = @(A,f,t,ph) A*sin(2*pi*f*t + ph);
dur = 2;
refRate = 120;
tStim = -0.5:1/refRate:dur;
stimAmp = 1;
stimFreq = 20;
stimFreq2 = [];
phase = 0;
phase2 = 0;
Fs = 1000;
tResponse = -0.5:1/Fs:dur;
noiseStd = 0;
plotFigs = 1;
%% STIMULUS
%% constructed neural response - use 120 Hz
% assumes that neurons respond to both onsets and offsets
% so if a stimulus at 60 Hz is [1 0 0], the neural response would be
% [1 1 0]. sampled at 120 Hz, this would be [1 0 1 0 0 0].
% 30 Hz
% stim = repmat([1 0 0 0],1,301/4);
% 40 Hz
% stim = repmat([1 0 1 0 0 0],1,301/6);
% stim(301) = 1;
%% make stimulus time series
stim = sinewave(stimAmp,stimFreq,tStim,phase);
if ~isempty(stimFreq2)
stim = stim + sinewave(stimAmp,stimFreq2,tStim,phase2);
end
stim(tStim<0) = 0; % zero baseline
if plotFigs % && 0
f(1) = figure;
plot(tStim, stim, '.-')
xlabel('time (s)')
ylabel('amp')
title(sprintf('%.2f Hz', stimFreq))
end
%% wavelet on stim
width = 10;
[spectrum,freqoi,timeoi] = ft_specest_wavelet(stim, tStim, 'freqoi', stimFreq, 'width', width);
specAmp = abs(squeeze(spectrum))';
if plotFigs
figure
plot(tStim, specAmp)
xlabel('time (ms)')
ylabel('spec amp')
title('wavelet')
end
% foi = 1:50;
% [spectrum,freqoi,timeoi] = ft_specest_wavelet(stim, tStim, 'freqoi', foi);
% specAmp = abs(squeeze(spectrum));
%
% if plotFigs
% ytick = 10:10:numel(freqoi);
% xtick = 51:50:numel(timeoi);
% figure
% imagesc(specAmp)
% rd_timeFreqPlotLabels(timeoi,freqoi,xtick,ytick);
% end
%% mtmconvol on stim
taper = 'hanning'; % 'dpss'
toi = tStim;
foi = stimFreq;
t_ftimwin = 12 ./ foi;
tapsmofrq = foi * 0.8; % for dpss
tfAmps = [];
[spectrum,ntaper,freqoi,timeoi] = ft_specest_mtmconvol(stim, tStim, ...
'timeoi', toi, 'freqoi', foi, 'timwin', t_ftimwin, ...
'taper', taper, 'tapsmofrq', tapsmofrq, 'dimord', 'chan_time_freqtap');
specAmp = abs(squeeze(spectrum))';
if plotFigs
figure
plot(tStim, specAmp)
xlabel('time (ms)')
ylabel('spec amp')
title('mtmconvol')
end
%% time-frequency on stim
taper = 'hanning';
foi = 1:50;
t_ftimwin = 10 ./ foi;
toi = tStim;
tfAmps = [];
[spectrum,ntaper,freqoi,timeoi] = ft_specest_mtmconvol(stim, tStim, ...
'timeoi', toi, 'freqoi', foi, 'timwin', t_ftimwin, ...
'taper', taper, 'dimord', 'chan_time_freqtap');
specAmp = abs(squeeze(spectrum))';
if plotFigs
ytick = 10:10:numel(foi);
xtick = 1:30:numel(toi);
figure
imagesc(specAmp)
rd_timeFreqPlotLabels(toi,foi,xtick,ytick);
title('mtmconvol')
end
%% SIMULATED NEURAL RESPONSE
%% make neural response
% by resampling the stimulus -- right thing to do?
ssvef = resample(stim, Fs, refRate);
ssvef = ssvef(1:length(tResponse));
noise = noiseStd.*randn(size(tResponse));
response = ssvef + noise;
if plotFigs
figure(f(1));
hold on
plot(tResponse, ssvef, 'r')
end
%% wavelet on simulated neural response
[spectrum,freqoi,timeoi] = ft_specest_wavelet(response, tResponse, 'freqoi', stimFreq);
specAmp = abs(squeeze(spectrum))';
if plotFigs
figure
plot(tResponse, specAmp)
xlabel('time (ms)')
ylabel('spec amp')
end
%% time-frequency on simulated neural response
taper = 'hanning';
foi = 1:80;
t_ftimwin = 10 ./ foi;
toi = tResponse;
tfAmps = [];
[spectrum,ntaper,freqoi,timeoi] = ft_specest_mtmconvol(response, tResponse, ...
'timeoi', toi, 'freqoi', foi, 'timwin', t_ftimwin, ...
'taper', taper, 'dimord', 'chan_time_freqtap');
specAmp = abs(squeeze(spectrum))';
if plotFigs
ytick = 10:10:numel(foi);
xtick = 1:500:numel(toi);
figure
imagesc(specAmp)
rd_timeFreqPlotLabels(toi,foi,xtick,ytick);
end