function [a, b, F, avgforce, avgadc] = force_calibration_automatic(forcetrace_filename, hires_filename, trace_index, show_plots)
% Calibrate mapping of high-resolution readout from minitweezers to force data.
% 
% [a, b] = force_calibration_automatic(forcetrace_filename, hires_filename, trace_index, show_plots)
%
% Computes scaling factor 'a' and additive shift 'b' to map high-resolution integer ADC voltage measurements 'adc' to force data 'force' as
%
% force = a * adc + b
%
% where 'force' is in units of pN, and 'adc' is an unsigned integer voltage measurement ADC reading.
%
% ARGUMENTS
%
%   forcetrace_filename (string) - name of minitweezers file containing forces (e.g. 'forcedata/p5ab_001/a8A.txt')
%   hires_filename (string) - name of .mat file containing high-resolution fiber data (e.g. 'Processed_p5ab/fiber1.mat')
%   trace_index (integer) - index of trace within high-resolution fiber data file
%   show_plots (boolean) - if true, then plots will be generated showing alignment of high-resolution data to force data
%
% RETURNS
%   a (float) - scaling factor for high-resolution data
%   b (float) - additive shift for high-resolution data
%   F (Nx1 float array) - calibrated 50 kHz force trace
%
% NOTES
%
% Operation roughly follows this scheme:
%
% * Identifies the subset of minitweezers force data where recording was switched to 1 kHz or, if no switch is detected, prints out a warning message.
%   The force data is plotted in the top plot, and 10 Hz data (not used here) is shaded in gray.
% * Subsamples (averages) the 50 kHz ADC data down to 1 kHz.  If plotting is enabled, the 50 kHz data is the second plot from the top, and the downsampled data is third from top.
% * The downsampled ADC data and the minitweezers 1 kHz force data are remapped to the interval [0,1], and dropped samples in the force data are added back in as zeros.  
%   The downsampled ADC data is padded at the end with its mean value.  
% * These transformed signals are fed to the fast image alignment algorithm published here:
%
%   http://www.mathworks.com/matlabcentral/fileexchange/18401
%
%   The reference for this method is:
%
%   Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup,
%   "Efficient subpixel image registration algorithms," Opt. Lett. 33,
%   156-158 (2008).

%   This algorithm is implemented in 'dftregistration.m', which I obtained from the Mathworks URL above.
%
%   * The downsampled, decimated 50 kHz data is then rescaled and shifted to match the mean and standard deviation of the minitweezers force data, and plotted (black) aligned in time against the force data (red) in the bottom plot.  
%   The middle two plots also have the 50 kHz data shifted in time to match up with where the image alignment algorithm thinks it should fit.

% Load minitweezers force trace.
fprintf('loading force data from "%s"...\n', forcetrace_filename);
[cycle, A_PsdY, A_PsdSum, A_Iris, B_PsdY, B_PsdSum, B_Iris, X_force, Y_force, Z_force, Tension, A_dist_Y, B_dist_Y, time, status] = textread(forcetrace_filename, '%d %f %f %f %f %f %f %f %f %f %f %f %f %f %d', 'headerlines', 2, 'commentstyle', 'shell');

% Determine recording time interval for minitweezers force data.
dt = (time(end) - time(1)) / (cycle(end) - cycle(1));
fprintf('minitweezers cycle dt = %.5f s\n', dt);

% Load high-resolution ADC recording trace.
fprintf('loading high-resolution data from "%s"...\n', hires_filename);
catalogue = whos('-file', hires_filename);
fibername = catalogue(1).name;
data = load(hires_filename, fibername);
fiber = getfield(data, fibername);
clear data;
trace = fiber{trace_index};

% Find the first time in the force trace where we switch from 10 Hz to 1 kHz recording.
indices = find(cycle(2:end) - cycle(1:end-1) == 1);
if (length(indices) == 0)
  % All force data is at 10 Hz -- we can't use this.
  fprintf('All force data is at 10 Hz\n');
  return;
else
  % Identify region of force data sampled at 1 kHz.
  index1 = indices(1);
  index2 = indices(end);
  clear indices;
  fprintf('1 kHz recording detected between index %d (%.3f s) and cycle %d (%.3f s)\n', cycle(index1), time(index1), cycle(index2), time(index2));

  % Set downsampling ratio.
  downsample_ratio = 50;
  fprintf('Downsampling (by averaging) %d : 1\n', downsample_ratio);
end

% Downsample the 50 kHz to 1 kHz trace by averaging.
fprintf('Downsampling 50 kHz signal...\n');
L = floor(length(trace) / downsample_ratio);
downsampled_trace = 0.0 * trace(1:L);
for n = 1:downsample_ratio
  x = trace(n:downsample_ratio:end);
  downsampled_trace = downsampled_trace + x(1:L);
end
downsampled_trace = downsampled_trace / double(downsample_ratio);

%% Find registration between the two using subpixel registration algorithm.

% Pad downsampled trace with zeros.
X = Y_force(index1:index2); % calibrated force measurement at 1 kHz
Y = downsampled_trace;      % uncalibrated ADC measurement downsampled to 1 kHz

% Remap traces to be alinged span interval [0,1].
X = (X - min(X)) / (max(X) - min(X));
Y = (Y - min(Y)) / (max(Y) - min(Y));

% Fill missing data in X (remapped calibrated force measurement) with zeros.
Z = zeros(cycle(index2) - cycle(index1) + 1, 1); % temporary
Z(cycle(index1:index2) - cycle(index1) + 1) = X;
X = Z;
clear Z;

% Pad Y with zeros.
Y = [Y; mean(Y) * ones(length(X) - length(Y), 1)];

% Find optimal registration using image alignment algorithm.
result = dftregistration(fft(X), fft(Y), 0.1);
error = result(1);
diffphase = result(2);
net_row_shift = result(3);
net_col_shift = result(4);

% Convert net row shift into a start time.
time1 = net_row_shift * dt + time(index1);
time2 = (length(downsampled_trace) - 1) * dt + time1;

% Compute indices of original force trace that match aligned region.
indices = find(time >= time1 & time <= time2);

% Compute parameters for transformation: force = a x + b
a = std(Y_force(indices)) / std(downsampled_trace);
b = - a * mean(downsampled_trace) + mean(Y_force(indices));

% DEBUG
avgforce = mean(Y_force(indices));
avgadc = mean(downsampled_trace);

% Compute calibrated 50 kHz force trace.
F = a * trace + b;

%
% Generate plot of both traces if plotting is requested.
%

if show_plots
  clf;

  % Plot 1 kHz trace

  nx = 1;
  ny = 4;

  subplot(ny,nx,1);
  hold on;
  plot(cycle(:)*dt + time(1), Y_force(:), 'r.');
  % Set axis
  ax = [time(1) time(end) min(Y_force) max(Y_force)];
  axis(ax);
  % Shade out 10 Hz regions if we have 1 kHz data present.
  fill([time(1) time(1) time(index1) time(index1)], [ax(3) ax(4) ax(4) ax(3)], 0.5 * [1 1 1], 'linestyle', 'none');
  fill([time(index2) time(index2) time(end) time(end)], [ax(3) ax(4) ax(4) ax(3)], 0.5 * [1 1 1], 'linestyle', 'none');
  title(sprintf('minitweezers force data from %s', forcetrace_filename), 'interpreter', 'none');
  xlabel('time / s');
  ylabel('force / pN');

  % Plot 50 kHz trace

  subplot(ny,nx,2);
  markersize = 1;
  h = plot(linspace(time1, time2, length(trace)), trace, 'k.', 'markersize', markersize);
  % Set axis
  ax = [time(1) time(end) min(trace) max(trace)];
  axis(ax);
  title(sprintf('50 kHz trace from %s trace %d', hires_filename, trace_index), 'interpreter', 'none');
  set(gca, 'xtick', []);
  set(gca, 'ytick', []);

  % Plot downsampled 50 kHz trace

  subplot(ny,nx,3);
  plot(linspace(time1, time2, length(downsampled_trace)), downsampled_trace, 'k.');
  ax = [time(1) time(end) min(trace) max(trace)];
  axis(ax);
  title(sprintf('50 kHz downsampled (by averaging) by ratio %d : 1', downsample_ratio));
  set(gca, 'xtick', []);
  set(gca, 'ytick', []);

  % Plot downsampled 50 kHz trace adjusted to axes

  subplot(ny,nx,4);
  hold on;
  plot(cycle(:)*dt + time(1), Y_force(:), 'r.');
  plot(linspace(time1, time2, length(downsampled_trace)), a*downsampled_trace + b, 'k.');
  ax = [time(1) time(end) min(Y_force) max(Y_force)];
  axis(ax);
  title(sprintf('50 kHz downsampled (by averaging) by ratio %d', downsample_ratio));
  xlabel('time / s');
  ylabel('force / pN');

  % Write transformation fit
  fprintf('transformation: force = %f x + %f\n', a, b);
end

return

