function [x,f,T] = compute_fourierCoefficients_matlab(s, win, noverlap, f, fs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculates a single complex fourier coefficient at frequency f and 
% returns the complex value
%
%   input:  s: time domain signal
%           f: frequency of fourier coefficients, in Hz
%           fs: sampling frequency of signal s
%           win: vector containing window function
%           noverlap: overlap in samples
%           
%
%   output: x: complex fourier coefficients
%           f: frequencies of fourier coefficients, same as input f
%           T: time in sec of window positions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 5
    fs = 1;
end

if nargin < 4
    error('Missing input arguments!');
end

s = s(:);
f = f(:);

win_len = length(win);
hopsize = win_len - noverlap;

t = ((0:1:(win_len-1))./fs)';
% t = ((-floor(win_len/2):1:floor(win_len/2))./fs)';
win_num = fix((length(s)-noverlap)/(win_len-noverlap));
x = zeros(win_num,length(f));
T = (win_len/2 : hopsize : length(s)-win_len/2)/fs;

twoPiT = 2*pi*t;

for f0 = 1:length(f)
%     fprintf(' Computing coefficient at %4.2f Hz.\n', f(f0));

    twoPiFt = f(f0)*twoPiT;
    cosine = cos(twoPiFt);
    sine = sin(twoPiFt);

    for w = 1:win_num
        start = (w-1)*hopsize+1;
        stop =  start + win_len -1;
        
        sig = s(start:stop) .*win;
        co = sum(sig .* cosine);
        si = sum(sig .* sine);
        x(w,f0) = (co + i*si);
    end
end
x = x';


%% compare
% a = angle(spectrogram(s,win,noverlap,f,fs));
% figure; imagesc(abs(angle(x)-a))