%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  Lock In Amp %%%%%% based on Hadas Benisty's code
%%%%%%%%%%%%%%%%%%%%%%%%%
function [Amp,theta,X,Y]=lockin(fs,f2,phiphi,x)
%% 
% fs - sampling frequency
% f2 - lockin frequency
% x  - the data could be a matrix (broken down by columns)
% phiphi - phase
phi=phiphi*pi/180; % is in rads; input is in degrees
%% start work

Nsamps = length(x(:,1)); % Number of samples
t = (0:Nsamps-1)'/fs; % Time vector
% correlator 
xc = (cos(f2*2*pi*t+phi)+1j*sin(2*pi*f2*t+phi))*2/Nsamps;
% multiply (integrate)
Cval = xc'*x; % Make sure data matrix is what it needs to be
%% extract amp and phase
Amp = abs(Cval); % Amplitude
theta = angle(Cval)*180/pi; % Phase also in degrees
X=real(Cval); % Inphase channel
Y=imag(Cval); % Quadrature channel