Python numpy.unwrap() Examples

def smoothness_score(angles):
    """Computes the smoothness score of a line interpolation according to the
    angles of each line.

    Args:
        - angles: Array of shape (n_interpolations, n_lines_per_interpolation)
            giving the angle of each line in each interpolation.

    Returns:
        - smoothness_scores: Array of shape (n_interpolations,) giving the
            average smoothness score for all of the provided interpolations.
    """
    angles = np.atleast_2d(angles)
    # Remove discontinuities larger than np.pi
    angles = np.unwrap(angles)
    diffs = np.abs(np.diff(angles, axis=-1))
    # Compute the angle difference from the first and last point
    total_diff = np.abs(angles[:, :1] - angles[:, -1:])
    # When total_diff is zero, there's no way to compute this score
    zero_diff = (total_diff < 1e-4).flatten()
    normalized_diffs = diffs/total_diff
    deviation = np.max(normalized_diffs, axis=-1) - 1./(angles.shape[1] - 1)
    # Set score to NaN when we aren't able to compute it
    deviation[zero_diff] = np.nan
    return deviation 

def phase_from_frequencyseries(htilde, remove_start_phase=True):
    """Returns the phase from the given frequency-domain waveform. This assumes
    that the waveform has been sampled finely enough that the phase cannot
    change by more than pi radians between each step.

    Parameters
    ----------
    htilde : FrequencySeries
        The waveform to get the phase for; must be a complex frequency series.
    remove_start_phase : {True, bool}
        Subtract the initial phase before returning.

    Returns
    -------
    FrequencySeries
        The phase of the waveform as a function of frequency.
    """
    p = numpy.unwrap(numpy.angle(htilde.data)).astype(
            real_same_precision_as(htilde))
    if remove_start_phase:
        p += -p[0]
    return FrequencySeries(p, delta_f=htilde.delta_f, epoch=htilde.epoch,
        copy=False) 

def test_specgram_angle_phase_equivalent(self):
        freqs = self.freqs_specgram
        speca, freqspeca, ta = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='angle')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspeca, freqspecp)
        assert_array_equal(ta, tp)
        assert_allclose(np.unwrap(speca, axis=0), specp,
                        atol=1e-06) 

def test_specgram_complex_phase_equivalent(self):
        freqs = self.freqs_specgram
        specc, freqspecc, tc = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='complex')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspecc, freqspecp)
        assert_array_equal(tc, tp)
        assert_allclose(np.unwrap(np.angle(specc), axis=0), specp,
                        atol=1e-06) 

def phaseanalysis(firstharmonic, displayplots=False):
    print('entering phaseanalysis')
    analytic_signal = hilbert(firstharmonic)
    amplitude_envelope = np.abs(analytic_signal)
    instantaneous_phase = np.angle(analytic_signal)
    if displayplots:
        print('making plots')
        fig = pl.figure()
        ax1 = fig.add_subplot(311)
        ax1.set_title('Analytic signal')
        X = np.linspace(0.0, 1.0, num=len(firstharmonic))
        pl.plot(X, analytic_signal.real, 'k', X, analytic_signal.imag, 'r')
        ax2 = fig.add_subplot(312)
        ax2.set_title('Phase')
        pl.plot(X, instantaneous_phase, 'g')
        ax3 = fig.add_subplot(313)
        ax3.set_title('Amplitude')
        pl.plot(X, amplitude_envelope, 'b')
        pl.show()
        pl.savefig('phaseanalysistest.jpg')
    instantaneous_phase = np.unwrap(instantaneous_phase)
    return instantaneous_phase, amplitude_envelope 

def test_specgram_complex_phase_equivalent(self):
        freqs = self.freqs_specgram
        specc, freqspecc, tc = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='complex')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspecc, freqspecp)
        assert_array_equal(tc, tp)
        assert_allclose(np.unwrap(np.angle(specc), axis=0), specp,
                        atol=1e-06) 

def test_specgram_angle_phase_equivalent(self):
        freqs = self.freqs_specgram
        speca, freqspeca, ta = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='angle')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspeca, freqspecp)
        assert_array_equal(ta, tp)
        assert_allclose(np.unwrap(speca, axis=0), specp,
                        atol=1e-06) 

def act(self, obs):
        quad_to_obj_pos, quad_to_obj_rot = obs[:2]
        quad_to_obj_T = transformations.quaternion_matrix(np.r_[quad_to_obj_rot[3], quad_to_obj_rot[:3]])
        quad_to_obj_T[:3, 3] = quad_to_obj_pos
        obj_to_quad_T = transformations.inverse_matrix(quad_to_obj_T)

        if self.tightness == 1.0:
            des_offset_hra = self.target_hra
        else:
            offset = obj_to_quad_T[:3, 3]
            offset_hra = xyz_to_hra(offset)
            target_hra = self.target_hra.copy()
            offset_hra[-1], target_hra[-1] = np.unwrap([offset_hra[-1], target_hra[-1]])
            des_offset_hra = (1 - self.tightness) * offset_hra + self.tightness * target_hra
        des_offset = hra_to_xyz(des_offset_hra)
        des_obj_to_quad_T = transformations.rotation_matrix(des_offset_hra[2], np.array([0, 0, 1]))
        des_obj_to_quad_T[:3, 3] = des_offset
        self.pbvs_pol.target_to_obj_T = transformations.inverse_matrix(des_obj_to_quad_T)
        return self.pbvs_pol.act(obs) 

def test_specgram_angle_phase_equivalent(self):
        freqs = self.freqs_specgram
        speca, freqspeca, ta = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='angle')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspeca, freqspecp)
        assert_array_equal(ta, tp)
        assert_allclose(np.unwrap(speca, axis=0), specp,
                        atol=1e-06) 

def test_specgram_angle_phase_equivalent(self):
        freqs = self.freqs_specgram
        speca, freqspeca, ta = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='angle')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspeca, freqspecp)
        assert_array_equal(ta, tp)
        assert_allclose(np.unwrap(speca, axis=0), specp,
                        atol=1e-06) 

def real_spectrum(signal, axis=-1, **kwargs):

    try:
        import matplotlib.pyplot as plt
    except ImportError:
        import warnings

        warnings.warn("Matplotlib is required for plotting")
        return

    S = np.fft.rfft(signal, axis=axis)
    f = np.arange(S.shape[axis]) / float(2 * S.shape[axis])

    plt.subplot(2, 1, 1)
    P = dB(S)
    plt.plot(f, P, **kwargs)

    plt.subplot(2, 1, 2)
    phi = np.unwrap(np.angle(S))
    plt.plot(f, phi, **kwargs) 

def test_specgram_complex_phase_equivalent(self):
        freqs = self.freqs_specgram
        specc, freqspecc, tc = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='complex')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspecc, freqspecp)
        assert_array_equal(tc, tp)
        assert_allclose(np.unwrap(np.angle(specc), axis=0), specp,
                        atol=1e-06) 

def test_specgram_complex_phase_equivalent(self):
        freqs = self.freqs_specgram
        specc, freqspecc, tc = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='complex')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspecc, freqspecp)
        assert_array_equal(tc, tp)
        assert_allclose(np.unwrap(np.angle(specc), axis=0), specp,
                        atol=1e-06) 

def test_specgram_angle_phase_equivalent(self):
        freqs = self.freqs_specgram
        speca, freqspeca, ta = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='angle')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspeca, freqspecp)
        assert_array_equal(ta, tp)
        assert_allclose(np.unwrap(speca, axis=0), specp,
                        atol=1e-06) 

def unwrap(p, discont=numpy.pi, axis=-1):
    """Unwrap by changing deltas between values to 2*pi complement.

    Args:
        p (cupy.ndarray): Input array.
        discont (float): Maximum discontinuity between values, default is
            ``pi``.
        axis (int): Axis along which unwrap will operate, default is the last
            axis.
    Returns:
        cupy.ndarray: The result array.

    .. seealso:: :func:`numpy.unwrap`
    """

    p = cupy.asarray(p)
    nd = p.ndim
    dd = sumprod.diff(p, axis=axis)
    slice1 = [slice(None, None)]*nd     # full slices
    slice1[axis] = slice(1, None)
    slice1 = tuple(slice1)
    ph_correct = _unwrap_correct(dd, discont)
    up = cupy.array(p, copy=True, dtype='d')
    up[slice1] = p[slice1] + cupy.cumsum(ph_correct, axis=axis)
    return up 

def testInstantaneousFrequency(self, shape, axis):
    # Instantaneous Frequency in numpy
    phase_np = np.pi * (2 * np.random.rand(*shape) - 1)
    unwrapped_np = np.unwrap(phase_np, axis=axis)
    dphase_np = np.diff(unwrapped_np, axis=axis)
    # Append with initial phase
    s = [slice(None),] * unwrapped_np.ndim
    s[axis] = slice(0, 1)
    slice_np = unwrapped_np[s]
    dphase_np = np.concatenate([slice_np, dphase_np], axis=axis) / np.pi

    phase_tf = tf.convert_to_tensor(phase_np)
    with self.cached_session() as sess:
      dphase_tf = sess.run(spectral_ops.instantaneous_frequency(phase_tf,
                                                                time_axis=axis))
    self.assertAllClose(dphase_np, dphase_tf) 

def test_specgram_complex_phase_equivalent(self):
        freqs = self.freqs_specgram
        specc, freqspecc, tc = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='complex')
        specp, freqspecp, tp = mlab.specgram(x=self.y,
                                             NFFT=self.NFFT_specgram,
                                             Fs=self.Fs,
                                             noverlap=self.nover_specgram,
                                             pad_to=self.pad_to_specgram,
                                             sides=self.sides,
                                             mode='phase')

        assert_array_equal(freqspecc, freqspecp)
        assert_array_equal(tc, tp)
        assert_allclose(np.unwrap(np.angle(specc), axis=0), specp,
                        atol=1e-06) 

def getPoseInterpolator(self, poseTimes, poses, unwrap_rpy=True):
        if unwrap_rpy:
            poses = np.array(poses, copy=True)
            poses[:,3:6] = np.unwrap(poses[:,3:6],axis=0)

        if self.interpolationMethod in ['slinear', 'quadratic', 'cubic']:
            f = scipy.interpolate.interp1d(poseTimes, poses, axis=0, kind=self.interpolationMethod)
        elif self.interpolationMethod == 'pchip':
            f = scipy.interpolate.PchipInterpolator(poseTimes, poses, axis=0)
        return f 

def calc_inst_info(modes,samplerate):
    """
Calculate the instantaneous frequency, amplitude, and phase of
each mode.
    """

    amp=np.zeros(modes.shape,np.float32)
    phase=np.zeros(modes.shape,np.float32)
    f=np.zeros(modes.shape,np.float32)

    print("Mode 1:", len(modes), samplerate)

    for m in range(len(modes)):
        h=scipy.signal.hilbert(modes[m])
        print(len(modes[m]))
        print("Mean Amplitude of mode ", m, np.mean(np.abs(h)))
        print("Mean Phase of mode ", m, np.mean(np.angle(h)))
        phase[m,:]=np.angle(h)
        print("Frequ", np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[m])]]))/(2*np.pi)*samplerate)
        amp[m,:]=np.abs(h)
        phase[m,:]=np.angle(h)
        f[m,:] = np.r_[np.nan,
                      0.5*(np.angle(-h[2:]*np.conj(h[0:-2]))+np.pi)/(2*np.pi) * samplerate,
                      np.nan]
        print("Mean Frequ of mode ", m, np.mean(np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[0])]]))/(2*np.pi)*samplerate))

        #f(m,:) = [nan 0.5*(angle(-h(t+1).*conj(h(t-1)))+pi)/(2*pi) * sr nan];

    # calc the freqs (old way)
    #f=np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[0])]]))/(2*np.pi)*samplerate

    # clip the freqs so they don't go below zero
    #f = f.clip(0,f.max())

    return f,amp,phase 

def _angles(skeletons):
    dd = np.diff(skeletons,axis=1);
    angles = np.arctan2(dd[...,0], dd[...,1])
    
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        angles = np.unwrap(angles, axis=1);
    
    mean_angles = np.mean(angles, axis=1)
    angles -= mean_angles[:, None]
    
    return angles, mean_angles 

def compute_inputs(x, freqs, n_fft, n_frames, input_features, norm=False):
    if norm:
        norm_fn = instance_norm
    else:
        def norm_fn(x):
            return x
    freqs_tf = tf.constant(freqs, name="freqs", dtype='float32')
    inputs = {}
    with tf.variable_scope('real'):
        inputs['real'] = norm_fn(tf.reshape(
            tf.matmul(x, tf.cos(freqs_tf)), [1, 1, n_frames, n_fft // 2]))
    with tf.variable_scope('imag'):
        inputs['imag'] = norm_fn(tf.reshape(
            tf.matmul(x, tf.sin(freqs_tf)), [1, 1, n_frames, n_fft // 2]))
    with tf.variable_scope('mags'):
        inputs['mags'] = norm_fn(tf.reshape(
            tf.sqrt(
                tf.maximum(1e-15, inputs['real'] * inputs['real'] + inputs[
                    'imag'] * inputs['imag'])), [1, 1, n_frames, n_fft // 2]))
    with tf.variable_scope('phase'):
        inputs['phase'] = norm_fn(tf.atan2(inputs['imag'], inputs['real']))
    with tf.variable_scope('unwrapped'):
        inputs['unwrapped'] = tf.py_func(
            unwrap, [inputs['phase']], tf.float32)
    with tf.variable_scope('unwrapped_difference'):
        inputs['unwrapped_difference'] = (tf.slice(
                inputs['unwrapped'],
                [0, 0, 0, 1], [-1, -1, -1, n_fft // 2 - 1]) -
            tf.slice(
                inputs['unwrapped'],
                [0, 0, 0, 0], [-1, -1, -1, n_fft // 2 - 1]))
    if 'unwrapped_difference' in input_features:
        for k, v in input_features:
            if k is not 'unwrapped_difference':
                inputs[k] = tf.slice(
                        v, [0, 0, 0, 0], [-1, -1, -1, n_fft // 2 - 1])
    net = tf.concat([inputs[i] for i in input_features], 1)
    return inputs, net 

def nanunwrap(x):
    '''correct for phase change for a vector with nan values     '''
    x = x.astype(np.float)

    bad = np.isnan(x)
    x = fillfnan(x)
    x = fillbnan(x)
    x = np.unwrap(x)
    x[bad] = np.nan
    return x 

def phase_from_polarizations(h_plus, h_cross, remove_start_phase=True):
    """Return gravitational wave phase

    Return the gravitation-wave phase from the h_plus and h_cross
    polarizations of the waveform. The returned phase is always
    positive and increasing with an initial phase of 0.

    Parameters
    ----------
    h_plus : TimeSeries
        An PyCBC TmeSeries vector that contains the plus polarization of the
        gravitational waveform.
    h_cross : TimeSeries
        A PyCBC TmeSeries vector that contains the cross polarization of the
        gravitational waveform.

    Returns
    -------
    GWPhase : TimeSeries
        A TimeSeries containing the gravitational wave phase.

    Examples
    --------s
    >>> from pycbc.waveform import get_td_waveform, phase_from_polarizations
    >>> hp, hc = get_td_waveform(approximant="TaylorT4", mass1=10, mass2=10,
                         f_lower=30, delta_t=1.0/4096)
    >>> phase = phase_from_polarizations(hp, hc)

    """
    p = numpy.unwrap(numpy.arctan2(h_cross.data, h_plus.data)).astype(
        real_same_precision_as(h_plus))
    if remove_start_phase:
        p += -p[0]
    return TimeSeries(p, delta_t=h_plus.delta_t, epoch=h_plus.start_time,
        copy=False) 

def _guess_delay(self,f_data,z_data):
        phase2 = np.unwrap(np.angle(z_data))
        gradient, intercept, r_value, p_value, std_err = stats.linregress(f_data,phase2)
        return gradient*(-1.)/(np.pi*2.) 

def plot_3D(self, args):
        def _plot(ax, key, val):
            if key == "view":
                for i in range(self.overlays):
                    x_url, d_url = self.d["entry/views/" + val].attrs[u"xy_" + str(i)].split(":")
                    si_x = self.si_prefix(self.d[x_url])
                    si_d = self.si_prefix(self.d[d_url])
                    x_index = int(self.d[x_url].shape[0] / 2)
                    y_index = int(self.d[x_url].shape[1] / 2)
                    ax.plot(self.d[x_url][x_index,y_index,:] / si_x[0], self.d[d_url][x_index,y_index,:] / si_d[0])
                    ax.set_ylabel("%s (%s%s)" % (self.d[d_url].attrs["name"], si_d[2], self.d[d_url].attrs["unit"]))
                    ax.set_xlabel("%s (%s%s)" % (self.d[x_url].attrs["name"], si_d[2], self.d[x_url].attrs["unit"]))
            elif key == "ds":
                if "/" not in val:
                    val = "entry/data0/" + val
                si_x = self.si_prefix(self.d[self.d[val].attrs["x_ds_url"]])
                si_y = self.si_prefix(self.d[self.d[val].attrs["y_ds_url"]])
                z_index = int(self.d[val].shape[2] / 2)
                data = self.d[val][:, :, z_index].T
                if self.unwrap_phase and val.split("/")[-1] == "phase":
                    data.T[~np.isnan(data.T)] = np.unwrap(data.T[~np.isnan(data.T)])
                if self.remove_offset_x_avg:
                    data -= np.nanmean(data, axis=1, keepdims=True)
                if self.remove_offset_y_avg:
                    data -= np.nanmean(data, axis=0, keepdims=True)
                ax.pcolorfast(self.d[self.d[val].attrs["x_ds_url"]][:self.d[val].shape[0]] / si_x[0], self.d[self.d[val].attrs["y_ds_url"]][:self.d[val].shape[1]] / si_y[0], data)
                ax.set_xlabel("%s (%s%s)" % (self.d[self.d[val].attrs["x_ds_url"]].attrs["name"], si_x[2], self.d[self.d[val].attrs["x_ds_url"]].attrs["unit"]))
                ax.set_ylabel("%s (%s%s)" % (self.d[self.d[val].attrs["y_ds_url"]].attrs["name"], si_y[2], self.d[self.d[val].attrs["y_ds_url"]].attrs["unit"]))
        ax_iter = iter(self.ax)
        for key, vals in args.items():
            for val in vals:
                _plot(ax_iter.next(), key, val) 

def phase_grad(self,spline = False):
        '''
        Unwrap, (optionally) spline smooth, and differentiate phase data. 
        The freuqnecy derivative of transmission or reflection phase data that follows 
        a arctan function yields a Lorentzian distribution.
        '''
        self.data = np.unwrap(self.data)
        if spline: self.spline_smooth()
        self.data = np.gradient(self.data)
        self._save_opt_data()

    #\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./=\./= 

def get_delay(self,f_data,z_data,delay=None,ignoreslope=True,guess=True):
		'''
		ignoreslope option not used here
		retrieves the cable delay assuming the ideal resonance has a circular shape
		modifies the cable delay until the shape Im(S21) vs Re(S21) is circular
		see "do_calibration"
		'''
		maxval = np.max(np.absolute(z_data))
		z_data = z_data/maxval
		A1, A2, A3, A4, fr, Ql = self._fit_skewed_lorentzian(f_data,z_data)
		if self.df_error/fr > 0.0001 or self.dQl_error/Ql>0.1:
			#print("WARNING: Calibration using Lorentz fit failed, trying phase fit...")
			A1 = np.mean(np.absolute(z_data))
			A2 = 0.
			A3 = 0.
			A4 = 0.
			#fr = np.mean(f_data)
			f = splrep(f_data,np.unwrap(np.angle(z_data)),k=5,s=self.phasefitsmooth)
			fr = f_data[np.argmax(np.absolute(splev(f_data,f,der=1)))]
			Ql = 1e4
		if ignoreslope==True:
			A2 = 0.
		else:
			A2 = 0.
			print("WARNING: The ignoreslope option is ignored! Corrections to the baseline should be done manually prior to fitting.")
			print("see also: resonator_tools.calibration.fit_baseline_amp() etc. for help on fitting the baseline.")
			print("There is also an example ipython notebook for using this function.")
			print("However, make sure to understand the impact of the baseline (parasitic coupled resonances etc.) on your system.")
			#z_data = (np.absolute(z_data)-A2*(f_data-fr)) * np.exp(np.angle(z_data)*1j)  #usually not necessary
		if delay is None:
			if guess==True:
				delay = self._guess_delay(f_data,z_data)
			else:
				delay=0.
			delay = self._fit_delay(f_data,z_data,delay,maxiter=200)
		params = [A1, A2, A3, A4, fr, Ql]
		return delay, params 

def _guess_delay(self,f_data,z_data):
        phase2 = np.unwrap(np.angle(z_data))
        gradient, intercept, r_value, p_value, std_err = stats.linregress(f_data,phase2)
        return gradient*(-1.)/(np.pi*2.) 

def pha(self, deg=False, unwrap=True, lag=True):
        """Phase of the electromagnetic field.

        Parameters
        ----------
        deg : bool
            If True the returned phase is in degrees, else in radians.
            Default is False (radians).

        unwrap : bool
            If True the returned phase is unwrapped.
            Default is True (unwrapped).

        lag : bool
            If True the returned phase is lag, else lead defined.
            Default is True (lag defined).

        """
        # Get phase, lead or lag defined.
        if lag:
            pha = np.angle(self.view())
        else:
            pha = np.angle(np.conj(self.view()))

        # Unwrap if `unwrap`.
        # np.unwrap removes the EMArray class;
        # for consistency, we wrap it in EMArray again.
        if unwrap and self.size > 1:
            pha = EMArray(np.unwrap(pha))

        # Convert to degrees if `deg`.
        if deg:
            pha *= 180/np.pi

        return pha


# 2. Input parameter checks for modelling

# 2.a <Check>s (alphabetically) 

def test_unwrap(self):
        q = [0., 3690., -270., 690.] * u.deg
        out = np.unwrap(q)
        expected = (np.unwrap(q.to_value(u.rad)) * u.rad).to(q.unit)
        assert out.unit == expected.unit
        assert np.allclose(out, expected, atol=1*u.urad, rtol=0)
        with pytest.raises(u.UnitsError):
            np.unwrap([1., 2.]*u.m)
        with pytest.raises(u.UnitsError):
            np.unwrap(q, discont=1.*u.m)