function result = ellipseFit4HC(x,y,options)%ellipseFit4HC Estimates the ellipse
parameters, based on N pairs of% measured (noisy) data (x and y), together
with their statistical% uncertainties (standard errors).%% SYNTAX:% result
= ellipseFit4HC(x,y,options)%% INPUT:% x - N-dimensional vector of x-
signal measurements; % y - N-dimensional vector of y-signal
measurements;% options - structure with the following parameters:%
- verbose: logical variable, flag for TABLE output of the%
estimated parameters (default value is true);% - isplot: logical
variable, flag for graphical output% (default value is true);%
- coverageFactor: coverage factor for computing the inteval%
estimators (default value is 2);% - alpha: nominal significance
level for computing the% (1-alpha)-quantiles (coverage factors) of
the inteval% estimators (default value is alpha = [], i.e. the %
options.coverageFactor is used);% - tolerance:
tolerance for convergence (default is 1e-12);% - maxloops: maximum
number of iterations (default is 100);% - correlation: known value
of the correlation coefficient, % rho = corr(x,y), while assuming
var(x) = var(y) = sigma2,% (default value is rho = 0);%
- displconst: displacement constant (default value is %
633.3*1000/(4*pi));% - displunit: displacement constant (default
value is % picometers [pm]);% - smoothByFFT: smooth
the original measurements (x,y)% by using the significant FFT
frequencies.%% OUPUT:% result - structure with detailed information on
estimated% parameters.%% EXAMPLE 1: (Fit the ellipse for generated
measurements x and y)% alpha0 = 0; beta0 = 0; % true ellipse center
[0,0]% alpha1 = 1; beta1 = 0.75; % true amplitudes of x and y signals %
phi0 = pi/3; % phase shift% X = @(t) alpha0 + alpha1 *
cos(t);% Y = @(t) beta0 + beta1 * sin(t + phi0);% sigma = 0.05;
% true measurement error STD% N = 50; % No. of
observations (x,y)% Ncycles = 0.8; % No. of whole ellipse
cycles% phi = Ncycles*(2*pi)*sort(rand(N,1)); % true phases % x = X(phi)
+ sigma*randn(size(X(phi)));% y = Y(phi) + sigma*randn(size(Y(phi)));%
result = ellipseFit4HC(x,y);%% EXAMPLE 2: (Fit the ellipse for generated
measurements x and y)% alpha0 = 0; beta0 = 0; % true ellipse center
[0,0]% alpha1 = 1; beta1 = 1; % true amplitudes% phi0 = pi/10;
% phase shift% X = @(t) alpha0 + alpha1 * cos(t);% Y = @(t)
beta0 + beta1 * sin(t + phi0);% sigma = 0.05; N = 1000; phi = (2*pi)*
((1:N)')/N; % x = X(phi) + sigma*randn(size(X(phi)));% y = Y(phi) +
sigma*randn(size(Y(phi)));% result = ellipseFit4HC(x,y)%
disp(result.TABLE_Displacements)%% EXAMPLE 3: (Ellipse fit based on N = 100000
interferometer measurements)% load InterferometerData% x = data(:,1) /
max(data(:,1));% y = data(:,2) / max(data(:,1));% result =
ellipseFit4HC(x,y)% disp(result.TABLE_Displacements(1:20,:))%% REMARKS:%
In particular, ellipseFit4HC can be useful for uncertaity evaluation% of the
estimated phases and/or displacements, based on quadrature% homodyne
interferometer measurements (with the Heydemann Correction% applied). %%
The *Heydemann Correction* is used to evaluate the phase in homodyne%
interferometer applications to correct the interferometer% nonlinearities),
for more details see e.g. Koening et al (2014) and% Wu, Su and Peng (1996).%
% Here we assume that the measurement errors for x and y are% independent
(optionally correlated, with known correlation% coefficient rho), with zero
mean and common variance sigma^2. %% The standard deviation of the
measurement errors, sigma, is assumed% to be small, such that the
measurements are relatively close to the% true, however unobservable ellipse
curve, as is the case for typical% interferometry measurements.% %
Moreover, due to numerical stability of the algorithm, it is% reasonable to
consider normalized measurements (x,y), i.e. such that% the length of the
main semiaxis of the fitted ellipse is close to 1.% % The standard algebraic
parametrization of ellipse, (A,B,C,D,E,F), see% e.g. Chernov & Wijewickrema
(2013), is % A*x^2 + 2*B*x*y + C*y^2 + 2*D*x + 2*E*y + F = 0 % or
(xc,yc,a,b,alpha), in geometric parametrization of the ellipse, of% the
following form:% [(x-xc)*cos(alpha) + (y-yc)*sin(alpha)]^2 / a^2 +%
... [(x-xc)*sin(alpha) + (y-yc)*cos(alpha)]^2 / b^2 = 1. % where -pi/2 <
alpha < pi/2, xc, yc denote the coordinates of the% ellipse center, and a, b
are the ellipse semiaxis.% % Alternatively, one can define the ellipse also
in the following% parametric form:% x(t) = xc + a * cos(t) * cos(alpha)
- b * sin(t) * sin(alpha)% y(t) = yc + a * cos(t) * sin(alpha) + b * sin(t)
* cos(alpha).% % However, here we consider the following algebraic
parametrization of% the ellipse, (B,C,D,F,G), as it is typically used in the
field of% interferometry, see e.g. Wu, Su and Peng (1996):% x^2 + B*y^2
+ C*x*y + D*x + F*y + G = 0,% and the geometric parametrization,
(alpha_0,alpha_1,beta_0, beta_1,% phi_0) of the form:% x(phi) = alpha_0
+ alpha_1 * cos(phi)% y(phi) = beta_0 + beta_1 * sin(phi + phi_0).%
where -pi/2 < phi0 < pi/2 is the phase offset, alpha_0, beta_0% denote the
coordinates of the ellipse center (offsets), and alpha_1,% beta_1 are the
signal amplitudes.% % Fitting ellipse based on measurements (x,y) is a non-
linear problem.% The presented approach is based on an approximate method
which is% correct to the first order. It is an iterative procedure based on%
subsequent first order Taylor expansions (linearizations) of the%
originally nonlinear model.% % The algorithm estimates the locally
approximate BLUEs (Best Linear% Unbiased Estimators) of the ellipse
parameters (B,C,D,F,G), the BLUEs% of the true signal values mu and nu (the
values on the fitted% ellipse) of the measurements (x,y), together with their
covariance% matrix, as suggested in Koening et al (2014).% % This is
based on iterative linearizations of the originally nonlinear% model with
nonlinear constraints on the model parameters. For more% details see Kubacek
(1988, p.152).% % Based on that, ellipseFit4HC estimates also the geometric%
parameters (alpha_0, alpha_1, beta_0, beta_1, phi_0) and the% N-
dimensional vector of phases phi (and/or displacements) together% with their
standard uncertainties computed by using the delta method.%% REFERENCES: % [1]
Koening, R., Wimmer, G. and Witkovsky V.: Ellipse fitting by% linearized
nonlinear constrains to demodulate quadrature homodyne% interferometer
signals and to determine the statistical uncertainty% of the interferometric
phase. To appear in Measurement Science and% Technology, 2014.% [2] Kubacek,
L.: Foundations of Estimation Theory. Elsevier, 1988. % [3] Chien-Ming Wu, Ching-
Shen Su and Gwo-Sheng Peng. Correction of% nonlinearity in one-frequency
optical interferometry. % Measurement Science and Technology, 7 (1996), 520?
24. % [4] Chernov N., Wijewickrema S.: Algorithms for projecting points onto%
conics. Journal of Computational and Applied Mathematics 251 (2013) % 8?1. %
(c) Viktor Witkovsky (witkovsky@savba.sk)% Ver.: 31-Jul-2014 18:27:32%% CHECK THE
INPUTS AND OUTPUTStic;narginchk(2, 3);if nargin < 3, options = []; endif
~isfield(options, 'smoothDataByFFT') options.smoothDataByFFT = false;endif
~isfield(options, 'displconst') options.displconst = 633.3*1000/(4*pi); %
unit: [pm] %options.displconst = 633.3/(4*pi); % unit: [nm] enddisplconst =
options.displconst;if ~isfield(options, 'displunit') options.displunit =
'picometer [pm]'; %options.displunit = 'nanometer [nm]';enddisplunit =
options.displunit;if ~isfield(options, 'tolerance') options.tolerance = 1e-
12;endtol = options.tolerance;if ~isfield(options, 'alpha') options.alpha =
[];endif ~isfield(options, 'maxloops') options.maxloops = 100;endif
~isfield(options, 'correlation') options.correlation = 0;endif
~isfield(options, 'coverageFactor') options.coverageFactor = 2;endif
~isfield(options, 'verbose') options.verbose = true;endverbose =
options.verbose;if ~isfield(options, 'isplot') options.isplot = true;endisplot
= options.isplot;%% Initialization of the Linear model with Type II constrainsn
= length(x);x = x(:);y = y(:);x_data = x;y_data = y;% rho is
know correlation coefficient rho = corr(x,y) [default rho = 0]% and var(x) =