% Example for a circular-linear fit of a bivariate Gaussian
% (by Eric Reifenstein, 1-June-2016)

% generating data for position and theta phase from a bivariate Gaussian
mu = [20,pi]; % vector of means
sigma = [36,-5;-5,1.2]; % covariance matrix
r = mvnrnd(mu,sigma,100); % draw 100 random samples
pos = r(:,1); % position data, in cm
phase = mod(r(:,2),2*pi); % theta phase, restricted to [0,2*pi[

% slope and offset:
% objective function, see Kempter et al. (2012) for details
myfun1 = @(p) -sqrt((sum(cos(phase-(p(1)*pos)))/length(phase)).^2 + (sum(sin(phase-(p(1)*pos)))/length(phase)).^2);
% finding the optimal slope, note that we have to restrict the range of
% possible slopes (here: one theta cycle per field width)
sl = fminbnd(myfun1, -2*pi/(max(pos)-min(pos)), 2*pi/(max(pos)-min(pos)));
% calculate offset
offs = atan2(sum(sin(phase-sl*pos)),sum(cos(phase-sl*pos))) + 2*pi;

% circular-linear correlation:
pos_circ = mod(abs(sl)*pos, 2*pi); % circular variable derived from the position
phase_mean = mod(angle(sum(exp(i*phase))/length(phase)),2*pi); % circular mean of the theta phase
pos_circ_mean = mod(angle(sum(exp(i*pos_circ))/length(phase)),2*pi); % circular mean of the circular position variable
% calculating the correlation
cc = sum(sin(phase - phase_mean) .* sin(pos_circ - pos_circ_mean)) / sqrt( sum(sin(phase - phase_mean).^2) * sum(sin(pos_circ - pos_circ_mean).^2) );

% plot the result, phases in deg
x = 0:0.1:max(pos)+10;
figure
plot(pos,180/pi * phase,'+')
hold on
plot(x, 180/pi * (sl * x + offs) )
xlabel('Position (cm)')
ylabel('Phase (deg)')
axis([0 inf 0 360])
set(gca,'YTick',[0:45:360])

% text output
fprintf('------------------------------------------------- \n')
fprintf('Example for a circular-linear fit \n')
fprintf('of a bivariate Gaussian \n')
fprintf('(by Eric Reifenstein, 1-June-2016) \n\n')
fprintf(sprintf('Slope:                         sl = %.1f deg/cm \n', 180/pi*sl))
fprintf(sprintf('Offset at position x=0:      offs = %.1f deg \n', 180/pi*offs))
fprintf(sprintf('Circular-linear correlation:   cc = %.3f \n', cc))
fprintf('------------------------------------------------- \n')
