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