function [TFmag TFphase faxis] = loadme
%
%
%--------------------------------------------------------------------------
[filename, pathname, filterindex] = uigetfile('*', 'PREsynaptic membrane potential..');
if (filename == 0),  return; end;
fullfname = sprintf('%s%s', pathname, filename);
pre  = load(fullfname);
cwd = pwd;
cd(pathname);
[filename, pathname, filterindex] = uigetfile('*', 'POSTsynaptic membrane potential..');
cd(cwd);
if (filename == 0),  return; end;
post = load(sprintf('%s%s', pathname, filename));
%--------------------------------------------------------------------------
Lpre  = length(pre);    Lpost = length(post);
if (Lpre ~= Lpost),     f = warndlg('Lengths of waveforms do not coincide', 'Error!'); 
 return; 
end;
%--------------------------------------------------------------------------
name          = 'Fourier Domain Preprocessing';
prompt        = 'Enter the Sampling Rate [kHz]';
numlines      = 1;
defaultanswer = {'10'};
answer        = inputdlg(prompt,name,numlines,defaultanswer);
dt            = 1./str2num(answer{1}); % ms --> 10kHz
time          = 0:dt:((Lpre * dt) - dt);
%--------------------------------------------------------------------------
pre(1)  = pre(2);
post(1) = post(2);
figure(1);  
subplot(2,1,1); plot(time(2:end)./1000., 1000*pre(2:end), 'k');
ylabel('pre [mV]', 'FontName', 'Arial', 'FontSize', 20);
set(gca, 'FontName', 'Arial', 'FontSize', 15, 'Box', 'on', 'XTickLabel', '');
%title('Data loaded');
set(gca, 'FontName', 'Arial', 'FontSize', 15, 'Box', 'on');
subplot(2,1,2); plot(time(2:end)./1000., 1000*post(2:end), 'r');
ylabel('post [mV]', 'FontName', 'Arial', 'FontSize', 20);
xlabel('time [s]', 'FontName', 'Arial', 'FontSize', 20);
set(gca, 'FontName', 'Arial', 'FontSize', 15);

rev = questdlg('Do you want to reverse pre <--> post ?', 'Reverse', 'Yes', 'No', 'Yes');
if (strcmp(rev, 'Yes')), tmp = post;, post = pre; pre = tmp; end;
%--------------------------------------------------------------------------

TFmag   = [];       TFphase = [];       faxis   = [];
[TFmag TFphase faxis] = tfpreprocessing(0., 0., pre, post, time);   
      


%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------

function [TFmag, TFphase, faxis] = tfpreprocessing(Tdiscard, Tdiscardsup, in, out, t)
%
% PREPROCESSING
% [TFmag, TFphase, faxis] = tfpreprocessing(Tdiscard, Tdiscardsup, in, out, t)
%
% Returns the magnitude and the phase of the complex transfer function
% in the frequency domain, between the input 'in' and the output 'out'
%
% Tdiscard    : [ms] how much initial interval to discard
% Tdiscardsup : [ms] how much final interval to discard
% in          : input waveform
% out         : output waveform
% t           : [ms] time axis
%
% Michele Giugliano and Corrado Cali', 2006, EPFL - Lausanne.
%
%--------------------------------------------------------------------------
dt = (t(21) - t(20)) * 1e-3; % sampling interval in [s]..
ind_cut     = max(find(t <= Tdiscard));
ind_cut_sup = max(find(t <= (t(end)-Tdiscardsup)));
i_input     = in(ind_cut:ind_cut_sup);
o_output    = out(ind_cut:ind_cut_sup);
clear in;   % frees some memory..
clear out;  % frees some memory..
%clear t;    % frees some memory..
%--------------------------------------------------------------------------
i_input  = i_input    - mean(i_input);  % remove the DC offset..
o_output = o_output   - mean(o_output); % remove the DC offset..

Tstart   = t(ind_cut) * 1e-3;        % [s]
Tfinal   = t(ind_cut_sup) * 1e-3;    % [s]
time     = Tstart:dt:Tfinal';        % New time axis column..
N        = length(i_input);          % Total number of samples..
%h       = hamming(N);               % symmetric Hamming window (FFT)
h = 1.;
i_input = i_input.*h;               % windowing the input..
o_output= o_output.*h;              % windowing the output..
clear time;                         % frees some memory..
clear Tstart;                       % frees some memory..
clear Tfinal;                       % frees some memory..
%--------------------------------------------------------------------------
NN      = 2^nextpow2(N);            % Next power of 2, for the FFT..
omega   = 1./(NN*dt);               % Frequency sampling interval (Hz)
faxis   = omega*(0:NN-1)';          % Frequency axis column..
Yin     = fft(i_input);             % Evaluating the FFT of the input
Yout    = fft(o_output);            % Evaluating the FFT of the output
Yin(1)  = nan;                      % Remove again the DC component..
Yout(1) = nan;                      % Remove again the DC component..
M       = NN / 2;                   % The FFT has a symmetry, you know..
faxis   = faxis(1:(M+1));           % Construct the frequency axis..
%--------------------------------------------------------------------------
Mag_in  = abs(Yin(1:(M+1)));        % Evaluate the magnitude of the in
Mag_out = abs(Yout(1:(M+1)));       % Evaluate the magnitude of the out
%--------------------------------------------------------------------------
Phase_in  = angle(Yin(1:(M+1)));    % Evaluate the phase of the in [rad]
Phase_out = angle(Yout(1:(M+1)));   % Evaluate the phase of the out [rad]
Phase_in  = Phase_in * (360./(2*pi)); % Convert it to degrees
Phase_out = Phase_out* (360./(2*pi)); % Convert it to degrees
%--------------------------------------------------------------------------
TFmag     = Mag_out ./ Mag_in;      % By definition of trasfer funtion..
TFphase   = unwrap(mod(Phase_out-Phase_in, 360) - 360); %