-
Notifications
You must be signed in to change notification settings - Fork 0
/
rd_simSampledSinewave2.m
136 lines (116 loc) · 3.23 KB
/
rd_simSampledSinewave2.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
% rd_simSampledSinewave2.m
%% setup
sinewave = @(A,f,t,ph) A*sin(2*pi*f*t + ph);
dur = 2;
refRate = 60;
tStim = -0.5:1/refRate:dur;
stimAmp = 1;
stimFreq = 15;
stimFreq2 = 20;
phase = pi/2;
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
%% 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 assuming that neurons respond to differences
stimUp = resample(stim, Fs, refRate);
ssvef = abs(diff(stimUp));
% ssvef = resample(ssvef, Fs, refRate);
if length(ssvef) > length(tResponse)
ssvef = ssvef(1:length(tResponse));
elseif length(tResponse) > length(ssvef)
tResponse = tResponse(1:length(ssvef));
end
noise = noiseStd.*randn(size(tResponse));
response = ssvef + noise;
if plotFigs
scaleFactor = max(stim)./max(ssvef);
figure
hold on
plot(tStim, stim, '.-')
plot(tResponse, ssvef*scaleFactor, 'r')
xlabel('time (s)')
ylabel('amp')
title(sprintf('%.2f Hz', stimFreq))
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