|
@@ -0,0 +1,85 @@
|
|
|
+function [spec_d, spec_t, spec_f] = get_spectrogram( x, t, win_size, t_resolution )
|
|
|
+% Usage: [spec_d, spec_t, spec_f] = get_spectrogram( x, t, win_size, t_resolution )
|
|
|
+%
|
|
|
+% Calculating amplitude spectrogram from signal x and time-vector t.
|
|
|
+%
|
|
|
+% -- input form --
|
|
|
+% x: Raw EEG signal (1-D vector)
|
|
|
+% y: Time vector (in millisecond resolution)
|
|
|
+% win_size: Size of sliding moving window (default: 2^10)
|
|
|
+% t_resolution: Jump size of sliding moving window (unit: sec, default: 100 msec)
|
|
|
+%
|
|
|
+% 2019-09-10
|
|
|
+%
|
|
|
+
|
|
|
+%% Default condition
|
|
|
+if nargin < 5
|
|
|
+ t_resolution = .1;
|
|
|
+end; if nargin < 4
|
|
|
+ win_size = 2^10;
|
|
|
+end
|
|
|
+
|
|
|
+srate = round( 1/nanmean(diff(t)) );
|
|
|
+t_fft = [t(1)+(((win_size*.5)+1)/srate), t(end)-(((win_size*.5)+1)/srate)];
|
|
|
+t_vec = linspace( t_fft(1), t_fft(end), (diff(t_fft)/t_resolution) +1);
|
|
|
+
|
|
|
+%% Sliding window index bank
|
|
|
+idx_collection = single([]);
|
|
|
+for tIdx = [1, length(t_vec)]
|
|
|
+ idx_collection(tIdx,:) = ...
|
|
|
+ [max(find( t<(t_vec(tIdx)))) - win_size*.5,...
|
|
|
+ (max(find( t<(t_vec(tIdx)))) + win_size*.5-1) ];
|
|
|
+end
|
|
|
+idx_collection(:,1) = linspace( idx_collection(1,1), idx_collection(end,1), length(t_vec) );
|
|
|
+idx_collection(:,2) = linspace( idx_collection(1,2), idx_collection(end,2), length(t_vec) );
|
|
|
+idx_collection = round(idx_collection);
|
|
|
+
|
|
|
+short = find( ~[diff(idx_collection')' == (win_size-1)]);
|
|
|
+idx_collection(short,2) = idx_collection(short,2)+1;
|
|
|
+try hann = hanning( idx_collection(1,2)-idx_collection(1,1)+1 )';
|
|
|
+catch hann = ones([1, win_size]); end
|
|
|
+
|
|
|
+%%
|
|
|
+spec_d = [];
|
|
|
+spec_t = [];
|
|
|
+for tIdx = 1:length(t_vec)
|
|
|
+ % (1) Indexing
|
|
|
+% idx = max(find( t<(t_vec(tIdx)))) - win_size*.5 :...
|
|
|
+% max(find( t<(t_vec(tIdx)))) + win_size*.5 -1 ;
|
|
|
+% actual_t = t(idx);
|
|
|
+% if tIdx==1, han=hanning(length(idx))'; end
|
|
|
+
|
|
|
+ idx_list = idx_collection(tIdx,1):idx_collection(tIdx,2);
|
|
|
+ d = hann .* x(idx_list(1:length(hann))); % to prevent 1-2 point mismatch
|
|
|
+
|
|
|
+ [fft_d,spec_f]= fft_half(d, srate);
|
|
|
+ spec_d( size(spec_d,1)+1 , :) = fft_d;
|
|
|
+ spec_t( length(spec_t)+1 ) = mean(t(idx_list(1:length(hann))));
|
|
|
+end
|
|
|
+spec_d = abs(spec_d);
|
|
|
+
|
|
|
+end
|
|
|
+
|
|
|
+%% Subfunctions
|
|
|
+function [s,f]=fft_half(x,Fs)
|
|
|
+N=(length(x)); k=0:N-1; T=N/Fs; f=k/T;
|
|
|
+cutOff = ceil(N/2); f = f(1:cutOff);
|
|
|
+nTrials = size(x, 1);
|
|
|
+if nTrials == 1
|
|
|
+ s = fft(x)/N*2; % normalize the data
|
|
|
+ s = s(1:cutOff); % Single trial FFT
|
|
|
+else
|
|
|
+ s = [];
|
|
|
+ for trialIdx = 1:nTrials
|
|
|
+ s_temp = fft(x(trialIdx,:))/N*2;
|
|
|
+ s_temp = s_temp(1:cutOff); % Single trial FFT
|
|
|
+ s = [s; s_temp];
|
|
|
+ end
|
|
|
+end
|
|
|
+end
|
|
|
+function [idx,short] = hb_findIdx( range, fullData )
|
|
|
+short=false;
|
|
|
+idx = max(find( fullData < range(1)))+1:max(find(fullData<range(2)));
|
|
|
+if length(idx)==0, idx = 1:max(find(fullData<range(2))); end
|
|
|
+if range(2) > fullData(end), short=1; end
|
|
|
+end
|