Source code for deerlab.correctphase

import numpy as np
from scipy.optimize import fminbound



def correctphase(V, full_output=False, offset=False):
    r"""
    Phase correction of complex-valued data.

    Rotates the phase of complex-valued data ``V`` to minimize the imaginary component.
    Among the two phases that minimize the imaginary part, the one that gives a real
    part with a positive average is used.
    For two-dimensional datasets ``V2D``, e.g. from measurements with multiple scans,
    each slice ``V2D[:,i]`` is phase-rotated independently.
    If the ``offset`` parameter is ``True``, the function will correct for a potential
    non-zero mean imaginary offset. 

    Parameters
    ----------
    V : array_like, or list of array_like
        Complex-valued 1D or 2D signal.

    full_output : boolean, optional
        If ``True``, return additional output arguments. (default: ``False``)

    offset : boolean 
        Enables numerical phase correction while accounting for a non-zero mean imaginary offset. 
        By default disabled. 

    Returns
    -------
    Vr : ndarray
        Real part of the phase-corrected data.

    Vi : ndarray (only if ``full_output==True``)
        Imaginary part of the phase-corrected data.

    phase : float scalar or ndarray (only if ``full_output==True``)
        Fitted phase, or list of phases for 2D data, used for correction, in radians.

    """
    
    if not np.iscomplexobj(V):
        raise ValueError("Data set must be complex-valued.")

    data1d = V.ndim==1

    V_2d = V.copy()
    if data1d:
        V_2d = V_2d[:, np.newaxis]

    if offset:
        # Account for a potential imaginary component with non-zero mean
        def objfcn(phi):
            Vim_corr = np.imag(V_2d*np.exp(1j*phi))
            Vim_corr -= np.average(Vim_corr,axis=0)
            return np.sum(Vim_corr**2)
        # Find one of the minima numerically
        phimin = np.atleast_1d(fminbound(objfcn, 0, np.pi, xtol=1.74e-3))

    else:
        # The follwing determines the phase that minimizes the cost
        # function = sum of squares of imaginary part of V*exp(1j*phi)
        # This cost function has the analytical form
        #
        #   (A+B) + (B-A)*cos(2*phi) + C*sin(2*phi)
        #    = offset + amp*cos(2*phi-phi0)
        #
        # where
        #    A = sum_k real(V_k)^2 / 2
        #    B = sum_k imag(V_k)^2 / 2
        #    C = sum_k real(V_k)*imag(V_k)
        #
        #    offset = A+B
        #    amp = sqrt((B-A)^2+C^2)
        #    phi0 = atan2(C, B-A)
        #
        # The cost function has two minima:
        #    phi = phi0/2 + pi/2   and   phi = phi0/2 + 3*pi/2

        # Calculate phase that minimizes cost function
        Vr = np.real(V_2d)
        Vi = np.imag(V_2d)
        A = np.sum(Vr**2, axis=0)/2
        B = np.sum(Vi**2, axis=0)/2
        C = np.sum(Vr*Vi, axis=0)
        phi0 = np.arctan2(C, B-A)
        phimin = phi0/2 + np.pi/2  # one of the two minimizers
    
    # Apply phase rotation
    V_2d *= np.exp(1j*phimin)[None,:]
    
    # Pick minimizer that yields positive average of real part
    reAvg = np.average(V_2d.real, axis=0)
    idx = reAvg < 0
    phimin[idx] += np.pi
    V_2d[:,idx] = -V_2d[:,idx]

    # Assemble output
    if data1d:
        V_2d = np.squeeze(V_2d, axis=1)
    Vreal = np.real(V_2d)
    Vimag = np.imag(V_2d)

    if offset: 
        Vreal += np.abs(np.average(Vimag, axis=0))
        Vimag -= np.average(Vimag, axis=0)

    if full_output:
        return Vreal, Vimag, phimin
    else:
        return Vreal