Opticist.org

       The intensity of an interferogram can be expressed as 

$ I_{ij} = A_{ij} + B_{ij} cos (\phi_j + \delta_i)$ (1)

where I is the intensity of the interferogram, the subscript "$i$" denotes the ith phase shifted image ($i = 1,2,...,M$) and "$j$" denotes the individual pixel location in each image ($j = 1,2,...,M$).  In the equation, $A_{ij}$ is the background or mean intensity, $B_{ij}$ is the modulation amplitude, $\phi_j$ is the angular phase information, $\delta_i$ is the induced phase shift amount of each of $M$ ($M\geq3$) acquired frames, and $N$ is the total number of pixels in each image frame. 

       From Eq. (1), it is evident that there are ($2MN+N+M$) unknowns, i.e., $A_{ij}$, $B_{ij}$, $\phi_j$, and $\delta_i$, and the number of available equations is ($MN$).  Using a direct method to solve the unknowns is not possible; instead, an iterative method can be employed.  The iteration may have two or more steps in each cycle and different groups of unknowns can be released in relevant steps to make the problem solvable.  The basic concept of the iterative algorithm is performing an iteration process until the sum of the squared differences between the estimated theoretical intensities and the measured intensities converge to a pre-defined small threshold value.  A brief description of the algorithm, its relations with conventional algorithms and a detailed procedure of using the algorithm are presented below.  

     1. Universal algorithm for phase extractions from interferograms with known phase shifts

       The conventional phase shifting algorithms require phase shift amounts not only to be known but also to be uniformly spaced.  In these algorithms, it is assumed that the background intensities and the modulation amplitudes have only pixel-to-pixel (intra-frame variations), i.e., $A_{1j} = A_{2j} = ... = A_{Mj} = A_{j}$  and in Eq. (1).  With the known phase shift amount  i, the number of unknowns is reduced to (3$N$).  Since the number of available equations is still ($MN$), the system of equations is posed as determined or overdetermined if  .  In this case, the unknowns can be readily solved.

       When $\phi_i$ is known and uniformly spaced, an explicit formulation for calculation of phase distributions can be derived.  This is the typical principle of conventional phase shifting algorithms.  When the phase shifts are known but randomly spaced, however, an explicit formulation is generally unavailable.  In this case, the following classical least-squares method, a solver for determined and overdetermined system of equations, can be employed to solve the problem.

       Defining a new set of variables as $a_j = A_{ij}, b_j = B_{ij} cos\phi_j$ and $c_j = -B_{ij}sin\phi_j$ can simplify the intensity Eq. (1) as $I_{ij} = a_j + b_j cos \delta_i + c_j sin \delta_i$  (2)

       The least squares error between theoretical and experimental interferograms, $S_j$ , accumulated from all the images described by Eq. (2), can be written as 

  $ S_j = \sum_{i=1}^M(I^e_{ij} - I_{ij})^2 = \sum_{i=1}^M (a_j + b_j  cos \delta_i + c_j  sin \delta_i - I_{ij}^e)^2       (3)$

where $I^e_{ij}$ is the experimentally measured intensity of the interferogram.  The least-squares criteria required for three unknowns ($a_j, b_j,   and  c_j$) can be expressed as

$\frac{\partial S_j}{\partial a_j} = 0, \frac{\partial S_j}{\partial b_j} = 0, \frac{\partial S_j}{\partial c_j} = 0,     (4)$

which yields

$ \pmatrix{M& \sum^M_{i=1} cos\delta_i & \sum^M_{i=1} sin \delta_i\\ \sum^M_{i=1} cos \delta_i & \sum^M_{i=1} cos^2 \delta_i & \sum^M_{i=1} cos \delta_i sin \delta_i\\ \sum^M_{i=1} sin \delta_i & \sum^M_{i=1} sin \delta_i cos \delta_i & \sum^M_{i=1} sin^2 \delta_i} \pmatrix{ a_j \\ b_j \\ c_j} = \pmatrix {\sum^M_{i=1} I^e_{ij} \\ \sum^M_{i=1} I^e_{ij} cos \delta_i \\ \sum^M_{i=1} I^e_{ij} sin \delta_i}.      (5)$

       From Eq. (5), the unknowns $a_j, b_j and c_j$ can be solved.  Then, the phase $\phi_j$ can be determined as

  $ \phi_j = tan^{-1} (\frac{-c_j}{b_j}).     (6)$

       The above phase shifting algorithm can provide the least-squares determination of the phase distributions when the phase shift amounts are accurately known.  Unlike the conventional algorithms, the universal algorithm does not require that the phase shift amounts have to be evenly spaced; this is the origin of the terminology, $\textit{universal}$.  Moreover, the processing time of the universal algorithm is of the same order as that of the conventional algorithms.  It is noted that when $M = 4$ and $\delta_i = (i-1) \pi /2 (i = 1,2,3,4), the universal phase shifting algorithm yields the formulation identical to the most-widely-used conventional four-frame algorithm.

     2. Inverse algorithm for determinations of phase shifts from interferograms with known phase distributions 

       Using the universal phase shifting algorithm, phase distributions can be extracted if the phase shift amounts are known.  If phase distributions are known, the phase shift amounts can be determined in a similar but inverse way.  In the inverse algorithm, it is assumed that the background intensities and modulation amplitudes have frame-to-frame variations but they are constants within each individual frame (inter-frame variations), i.e., $A_{i1} = A_{i2} = ... = A{iN} = A_i$ and $B_{i1} = B_{i2} = ... = B_{iN} = B_i$  in Eq. (1).  Consequently, the number of unknowns in the equation is reduced to ($3M+N$).  If the phase   is known, the number of unknowns can be further reduced to (3$M$).  With the number of available equations still being ($MN$), a similar overdeterministic least-squares method can be employed to obtain the unknowns.

       Defining a new set of variables for each frame as $a'_i = A_{ij}, b'+i = B_{ij} cos \delta_i$ and $c'_i = -Bij sin delta_i$, the intensity Eq. (1) can be expressed as

$I_{ij} = a'_i + b'_i  cos  \phi_j + c'_i  sin  \phi_j.    (7)$

      The least squares error, $S'_i$ , accumulated from all the pixels in the $i$th image, can be expressed as

  $ S'_j  = \sum^N_{j=1} (I_{ij} - I^e_{ij})^2 = \sum^N_{j=1} (a'_i + b'_i  cos  \phi_j + c'_j  sin  \phi_j - I^e_{ij})^2.      (8)$

      When $\phi_j$ is known, the least squares criteria of the unknowns $a'_i, b'_j, and c'_j$ yield

$ \pmatrix{M& \sum^M_{i=1} cos\phi_i & \sum^M_{i=1} sin \phi_i\\ \sum^M_{i=1} cos \phi_i & \sum^M_{i=1} cos^2 \phi_i & \sum^M_{i=1} cos \phi_i sin \phi_i\\ \sum^M_{i=1} sin \phi_i & \sum^M_{i=1} sin \phi_i cos \phi_i & \sum^M_{i=1} sin^2 \phi_i} \pmatrix{ a_j \\ b_j \\ c_j} = \pmatrix {\sum^M_{i=1} I^e_{ij} \\ \sum^M_{i=1} I^e_{ij} cos \phi_i \\ \sum^M_{i=1} I^e_{ij} sin \phi_i}.         (9)$

      The unknowns $a'_i ,b'_i$ and $c'_i$ can be obtained from Eq. (9) and the phase shifts $\delta_i$ can be determined from the unknowns as

  $\delta_i = tan^-1(\frac{-c'_i}{b'_i}.     (10)$

      The above algorithm can be regarded as an inverse procedure of the universal phase shifting algorithm.  In a typical phase shifting analysis, the inverse algorithm can be employed to verify the phase shift amounts, and the new detected phase shift amounts can be used to adjust the input of phase shifts for phase shifting analysis.  This procedure can help increase the accuracies of phase extractions.  This is the concept of the AIA for random phase shifting.

    2.3 Iterative strategy for random phase shifting

      The combination of the two algorithms (universal phase shifting algorithm and inverse phase shifting algorithm) yields an iteration procedure that can determine phase distributions and phase shifts simultaneously from randomly phase-shifted interferograms.  The AIA contains three steps in each iteration cycle.  In the k-th iteration cycle:

      Step 1 - Employ the universal phase shifting algorithm, i.e., Eqs. (5) and (6), to calculate phase distributions $\phi^{k-1}_j(j = 1,2,...,N)$ based on phase shifts $\delta^{k - 1}_i (i = 1,2,...,M)$ obtained in the second step of previous iteration cycle, i.e., the (k-1)-th one.  For the first iteration cycle, an initial estimation of phase shifts is not required because they can be any random values. 

      Step 2 - With known phase distributions $\phi^{k-1}_j(j = 1,2,...,N)$ , utilize the inverse phase shifting algorithm, i.e., Eqs. (9) and (10), to determine updated phase shifts,  $\delta^{k}_i (i = 1,2,...,M)$.

      Step 3 - Check to see whether the iteration results satisfy the converging criteria.  It is the relative phase shift amounts that will converge, so the converging criteria are

$|(\delta^k_i  -  \delta^k_i)  -  (\delta^{k-1}_i  -  \delta^{k-1}_1)|  < \epsilon,   i = 2,..., M,            (11)$

where $\epsilon$ is the pre-defined threshold of accuracy, e.g., $10^{1-6}$.  If the converging criteria are satisfied, then $\phi^{k-1}_j(j = 1,2,...,N)$ and  $\delta^{k}_i (i = 1,2,...,M)$ are the final phase distributions and phase shifts, respectively.  Otherwise, the iteration will repeat until the criteria are satisfied.