def angle(x):
    # Wrapper for angle that handles CXData objects
    if isinstance(x, CXData):
        l=[]
        for i in xrange(len(x)):
            l.append(sp.angle(x.data[i]))

        return CXData(data=l)
    elif isinstance(x, np.ndarray):
        return sp.angle(x)
    else:
        raise Exception('Unknown data type passed to angle')

def symmetrydemo():
    a=sin(linspace(-5*pi,5*pi,10000))
    b=a+2
    c=a-0.5
    ah,bh,ch=hilbert(a),hilbert(b),hilbert(c)
    ph_a,ph_b,ph_c=unwrap(angle(ah)),unwrap(angle(bh)),unwrap(angle(ch))
    omega_a=diff(ph_a)
    omega_b=diff(ph_b)
    omega_c=diff(ph_c)
    subplot(211),plot(ph_a),plot(ph_b),plot(ph_c)
    subplot(212),plot(omega_a),plot(omega_b),plot(omega_c)
    grid()
    show()
    return a,b,c

 def get_fft(self, fs, taps, Npts):
     Ts = 1.0/fs
     fftpts = fftpack.fft(taps, Npts)
     self.freq = scipy.arange(0, fs, 1.0/(Npts*Ts))        
     self.fftdB = 20.0*scipy.log10(abs(fftpts))
     self.fftDeg = scipy.unwrap(scipy.angle(fftpts))
     self.groupDelay = -scipy.diff(self.fftDeg)

def ACphase(tf, unit='deg', unwrapTol=0.5):
	"""
	Return the phase in desired *unit* of a transfer function *tf*
	
	* ``deg`` stands for degrees
	* ``rad`` stands for radians
	
	The phase is unwrapped (discontinuities are stiched together to make it 
	continuous). The tolerance of the unwrapping (in radians) is 
	*unwrapTol* times ``pi``. 
	"""
	# Get argument
	ph=angle(tf)
	
	# Unwrap if requested
	if (unwrapTol>0) and (unwrapTol<1):
		ph=unwrap(ph, unwrapTol*pi)
	
	# Convert to requested unit
	if unit=='deg':
		return ph/pi*180.0
	elif unit=='rad':
		return ph
	else:
		raise Exception, "Bad phase unit."

def mmse_stsa(infile, outfile, noise_sum):
    signal, params = read_signal(infile, WINSIZE)
    nf = len(signal)/(WINSIZE/2) - 1
    sig_out=sp.zeros(len(signal),sp.float32)

    G = sp.ones(WINSIZE)
    prevGamma = G
    alpha = 0.98
    window = sp.hanning(WINSIZE)
    gamma15=spc.gamma(1.5)
    lambdaD = noise_sum / 5.0
    percentage = 0
    for no in xrange(nf):
        p = int(math.floor(1. * no / nf * 100))
        if (p > percentage):
            percentage = p
            print "{}%".format(p),

        y = get_frame(signal, WINSIZE, no)
        Y = sp.fft(y*window)
        Yr = sp.absolute(Y)
        Yp = sp.angle(Y)
        gamma = Yr**2/lambdaD
        xi = alpha * G**2 * prevGamma + (1-alpha)*sp.maximum(gamma-1, 0)
        prevGamma = gamma
        nu = gamma * xi / (1+xi)
        G = (gamma15 * sp.sqrt(nu) / gamma ) * sp.exp(-nu/2) * ((1+nu)*spc.i0(nu/2)+nu*spc.i1(nu/2))
        idx = sp.isnan(G) + sp.isinf(G)
        G[idx] = xi[idx] / (xi[idx] + 1)
        Yr = G * Yr
        Y = Yr * sp.exp(Yp*1j)
        y_o = sp.real(sp.ifft(Y))
        add_signal(sig_out, y_o, WINSIZE, no)
    
    write_signal(outfile, params, sig_out)

def plot_freqz(b, a, w = None, npoints = None, title = '', db = False, createFigure = True, label = ''):
    # Create the omega array if necessary
    if npoints is None:
        npoints = 1000

    if w is None:
        w = scipy.arange(-scipy.pi, scipy.pi, 2*scipy.pi/(npoints))

    # Calculate the frequency response
    d = loudia.freqz(b.T, a.T, w)

    if db:
        mag = 20.0 * scipy.log10(abs(d[:,0]))
    else:
        mag = abs(d[:,0])

    import pylab
    if createFigure:
        pylab.figure()
        
    pylab.subplot(2,1,1)
    pylab.plot(w, mag, label = label)
    pylab.title('%s \n Magnitude of the Frequency Response' % title)

    pylab.subplot(2,1,2)
    pylab.plot(w, scipy.angle(d[:,0]), label = label)
    pylab.title('Angle of the Frequency Response')

def csr3(complex_n):
    ang = sp.angle(complex_n) # sp.arctan(a.imag/a.real) why it does not work?!?!
    r = sp.sqrt(sp.square(complex_n.real)+sp.square(complex_n.imag))
    if (sp.sin(ang/2)>=0): #sin>0
        return sp.sqrt(r)*(complex(sp.cos(ang/2),sp.sin(ang/2)))
    else:
        return sp.sqrt(r)*(complex(sp.cos((ang/2)+sp.pi),sp.sin((ang/2)+sp.pi)))

def floquet_modes(H, T, args=None, sort=False):
    """
    Calculate the initial Floquet modes Phi_alpha(0) for a driven system with
    period T.
    
    Returns a list of :class:`qutip.Qobj` instances representing the Floquet
    modes and a list of corresponding quasienergies, sorted by increasing
    quasienergy in the interval [-pi/T, pi/T]. The optional parameter `sort`
    decides if the output is to be sorted in increasing quasienergies or not.

    Parameters
    ----------
    
    H : :class:`qutip.Qobj`
        system Hamiltonian, time-dependent with period `T`
        
    args : dictionary
        dictionary with variables required to evaluate H
     
    T : float
        The period of the time-dependence of the hamiltonian. The default value
        'None' indicates that the 'tlist' spans a single period of the driving.
     
    Returns
    -------

    output : list of kets, list of quasi energies

        Two lists: the Floquet modes as kets and the quasi energies.

    """

    # get the unitary propagator
    U = propagator(H, T, [], args)

    # find the eigenstates for the propagator
    evals,evecs = la.eig(U.full())

    eargs = angle(evals)
    
    # make sure that the phase is in the interval [-pi, pi], so that the
    # quasi energy is in the interval [-pi/T, pi/T] where T is the period of the
    # driving.
    #eargs  += (eargs <= -2*pi) * (2*pi) + (eargs > 0) * (-2*pi)
    eargs  += (eargs <= -pi) * (2*pi) + (eargs > pi) * (-2*pi)
    e_quasi = -eargs/T

    # sort by the quasi energy
    if sort == True:
        order = np.argsort(-e_quasi)
    else:
        order = list(range(len(evals)))

    # prepare a list of kets for the floquet states
    new_dims  = [U.dims[0], [1] * len(U.dims[0])]
    new_shape = [U.shape[0], 1]
    kets_order = [Qobj(np.matrix(evecs[:,o]).T, dims=new_dims, shape=new_shape) for o in order]

    return kets_order, e_quasi[order]

def getinstfreq(imfs):
    omega=zeros((imfs.shape[0],imfs.shape[1]-1),dtype=float)
    for i in range(imfs.shape[0]):
        h=hilbert(imfs[i,:])
        theta=unwrap(angle(h))
        omega[i,:]=diff(theta)
        
    return omega

 def phase(self, degrees=False):
     '''
     Return an array of the phase angles of the wavelet coefficients,
     in radians (set degrees=True for degrees).
     '''
     phase = sp.angle(self.wave)
     if degrees:
         phase *= 180 / sp.pi
     return phase

 def compute_by_noise_pow(self,signal,n_pow):
     s_spec = sp.fft(signal*self._window)
     s_amp = sp.absolute(s_spec)
     s_phase = sp.angle(s_spec)
     amp = s_amp**2.0 - n_pow*self._coefficient
     amp = sp.maximum(amp,0.0)
     amp = sp.sqrt(amp)
     amp = self._ratio*amp + (1.0-self._ratio)*s_amp
     spec = amp * sp.exp(s_phase*1j)
     return sp.real(sp.ifft(spec))

def steady_state_response(zs, rmin, rmax):
    """
    Returns a plot with the steady state response of a
    single degree of freedom damped system.

    Parameters
    ----------
    zs: array
        Array with the damping values
    rmin, rmax: float
        Minimum and maximum frequency ratio

    Returns
    ----------
    r: Array
        Array containing the values for the frequency ratio
    A: Array
        Array containing the values for anmplitude

        Plot with steady state magnitude and phase

    Examples:
    >>> r, A = steady_state_response([0.1, 0.3, 0.8], 0, 2)
    >>> A[10]
    (0.98423159842039087-0.15988334018879749j)
    """

    if not isinstance(zs, list):
        zs = [zs]
    r = sp.linspace(rmin, rmax, 100*(rmax-rmin))
    A0 = sp.zeros((len(zs), len(r)), complex)
    for z in enumerate(zs):
        A0[z[0]] = (1/(1 - r**2 + 2*1j*r*z[1]))

    fig = plt.figure()
    ax1 = fig.add_subplot(211)
    ax2 = fig.add_subplot(212, sharex=ax1)
    plt.tight_layout()

    ax1.set_ylabel('Normalized Amplitude (dB)')
    ax1.set_title('Normalized Amplitude vs Frequency Ratio')

    ax2.set_xlabel('Frequency Ratio')
    ax2.set_ylabel('Phase Lag (deg)')
    ax2.set_title('Phase vs Frequency Ratio')

    for A in A0:
        ax1.plot(r, (sp.absolute(A)))
        ax2.plot(r, -sp.angle(A)/sp.pi*180)

    ax1.legend((['$\zeta$ = ' + (str(s)) for s in zs]))

    _ = plt.show()

    return r, A

    def init_data(self, *args, **kwargs):

        if args[0] == 'det_mod':
            if CXP.actions.preprocess_data:
                self.read_in_data()
            else:
                self.load()

        elif args[0] == 'probe_det_mod':
            if CXP.actions.preprocess_data:
                #  Get list of white files
                CXP.log.info('Preprocessing probe detector modulus.')
                if CXP.io.whitefield_filename not in [None, '']: # If whitefields were measured

                    wfilename, wfilerange, wn_acqs = [CXP.io.whitefield_filename, CXP.io.whitefield_filename_range,
                                                      CXP.measurement.n_acqs_whitefield]
                    self.pattern = wfilename.count('{')
                    if self.pattern == 1:
                        wf = [wfilename.format(i) for i in range(wfilerange[0], wfilerange[1])]
                    elif self.pattern == 2:
                        wf = [wfilename.format(wfilerange[0], i) for i in range(wn_acqs)]
                    elif self.pattern == 3:
                        wf = glob.glob(wfilename.split('}')[0]+'}*')

                    res = self.preprocess_data_stack(0, 1, wf, self.pattern, None, None, no_decorate=True)
                    self.data = res[1]
                else: #Guesstimate the whitefield from the average of the diffraction patterns
                    pass
            else:
                self.load(CXData.raw_data_filename_string.format('probe_det_mod'))

            try:
                probe = self.__class__.__all__['probe']
                probe.data[0] = spf.ifft2(self.data[0]*exp(complex(0., 1.)*sp.angle(spf.fft2(probe.data[0]))))
                CXP.log.info('Applied probe modulus constraint.')
            except (AttributeError, KeyError):
                pass

        elif args[0] == 'dark':
            if CXP.actions.preprocess_data:
                # Get list of dark files
                CXP.log.info('Preprocessing darkfield.')
                dfilename, dfilerange, dn_acqs = [CXP.io.darkfield_filename, CXP.io.darkfield_filename_range,
                                                  CXP.measurement.n_acqs_darkfield]
                self.pattern = dfilename.count('{')
                if self.pattern == 1:
                    df = [dfilename.format(i) for i in range(dfilerange[0], dfilerange[1])]
                elif self.pattern == 2:
                    df = [dfilename.format(dfilerange[0], i) for i in range(dn_acqs)]
                elif self.pattern == 3:
                    df = glob.glob(dfilename.split('}')[0]+'}*')
                res = self.preprocess_data_stack(0, 1, df, self.pattern, None, None, no_decorate=True)
                self.data = res[1]
            else:
                self.load(CXData.raw_data_filename_string.format('probe_det_mod'))

 def assert_is_analytic(self, signal, amlaw=None):
     """Assert that signal is analytic."""
     omega = angle(signal)
     if amlaw is not None:
         recons = np.exp(1j * omega) * amlaw
     else:
         recons = np.exp(1j * omega)
     real_identical = np.allclose(np.real(recons), np.real(signal))
     imag_identical = np.allclose(np.imag(recons), np.imag(signal))
     if not (imag_identical and real_identical):
         raise AssertionError("Signal is not analytic.")

 def test_anaqpsk(self):
     """Test quaternary PSK signal."""
     signal, phases = ana.anaqpsk(512, 64, 0.25)
     self.assert_is_analytic(signal)
     # Count discontinuities in the signal and the phases and assert that
     # they appear in the same locations
     uphase = unwrap(angle(signal))
     dphase = np.diff(uphase)
     base_value = mode(dphase)[0][0]
     signal_phase_change = np.abs(dphase - base_value) > 0.0001
     ideal_phase_change = np.diff(phases) != 0
     np.testing.assert_allclose(signal_phase_change, ideal_phase_change)

def hilbert_phaser(s):
    """
    FUNC: hilbert_phaser
    DESCR: Return the instantaneous phase of s using the hilbert transform
    """
    #f = file("hilbert_fpe.pickle", "w")
    #pickle.dump(s, f)
    #f.close()
    #print "utils:hilbert_phaser(): sig.hilbert(", s[0:10], ")"
    h = scipy.signal.hilbert(s)
    #print "utils:hilbert_phaser(): scipy.absolute(", h[0:10], ")"    
    a = scipy.absolute(h)
    #print "utils:hilbert_phaser(): scipy.angle(", h[0:10], ")"    
    phase = scipy.angle(h)
    return phase

def calc_noise(filepath):
    noise_sum=None
    signal, params = read_signal(filepath, WINSIZE)
    nf = len(signal)/(WINSIZE/2) - 1
    noise_sum=sp.zeros(WINSIZE,sp.float32)
    window = sp.hanning(WINSIZE)
    for no in xrange(nf):
        y = get_frame(signal, WINSIZE, no)
        Y = sp.fft(y*window)
        Yr = sp.absolute(Y)
        Yp = sp.angle(Y)
        if ( no < 20 ):
            noise_sum = noise_sum + Yr**2
        else:
            break
    return noise_sum

def frfplot(f, H):
    plt.subplot(211)
    plt.plot(f, 20 * sp.log10(sp.absolute(sp.sum(H, axis=1))), '-')
    plt.grid('on')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('FRF (dB)')
    axlim = plt.axis()
    plt.axis(
        axlim + sp.array([0, 0, -0.1 * (axlim[3] - axlim[2]), 0.1 * (axlim[3] - axlim[2])]))

    plt.subplot(212)
    plt.plot(f, sp.unwrap(sp.angle(sp.sum(H, axis=1))) / sp.pi * 180, '-')
    plt.grid('on')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Phase (deg)')
    axlim = plt.axis()
    plt.axis(
        axlim + sp.array([0, 0, -0.1 * (axlim[3] - axlim[2]), 0.1 * (axlim[3] - axlim[2])]))

def inst_freq(x, t=None, L=1):
    """
    Compute the instantaneous frequency of an analytic signal at specific
    time instants using the trapezoidal integration rule.

    :param x: The input analytic signal
    :param t: The time instants at which to calculate the instantaneous frequencies.
    :param L: Non default values are currently not supported.
        If L is 1, the normalized instantaneous frequency is computed. If L > 1,
        the maximum likelihood estimate of the instantaneous frequency of the
        deterministic part of the signal.
    :type x: numpy.ndarray
    :type t: numpy.ndarray
    :type L: int
    :return: instantaneous frequencies of the input signal.
    :rtype: numpy.ndarray
    :Example:
    >>> from tftb.generators import fmsin
    >>> x = fmsin(70, 0.05, 0.35, 25)[0]
    >>> instf, timestamps = inst_freq(x)
    >>> plot(timestamps, instf) #doctest: +SKIP

    .. plot:: docstring_plots/processing/freq_domain/inst_freq.py
    """
    if x.ndim != 1:
        if 1 not in x.shape:
            raise TypeError("Input should be a one dimensional array.")
        else:
            x = x.ravel()
    if t is not None:
        if t.ndim != 1:
            if 1 not in t.shape:
                raise TypeError("Time instants should be a one dimensional "
                                "array.")
            else:
                t = t.ravel()
    else:
        t = np.arange(2, len(x))

    fnorm = 0.5 * (angle(-x[t] * np.conj(x[t - 2])) + np.pi) / (2 * np.pi)
    return fnorm, t

        try :
            if (x.shape[0]==0) or (y.shape[0]== 0) : raise ValueError("Empty signal")
            if x.shape[0]!=y.shape[0] : raise ValueError("The two signals have different length")
        except ValueError, err_msg:
            raise ValueError(err_msg)
            return
        
    
        M=x.shape[0]
                

        #computing the analytic signal and the instantaneous phase
        x_analytic=hilbert(np.hstack(x.values))
        y_analytic=hilbert(np.hstack(y.values))
        
        phx=np.unwrap(scipy.angle(x_analytic))
        phy=np.unwrap(scipy.angle(y_analytic))
    
         
        disc_perc=np.floor(x.shape[0]/10.0)    
        
        phx_s=phx[disc_perc-1:M-disc_perc]
        phy_s=phy[disc_perc-1:M-disc_perc]
        
        ph_nm=(self._n*phx_s-self._m*phy_s)
        psi_nm=np.mod(ph_nm, 2*np.pi)

        gamma2_nm=np.square(np.mean(np.cos(psi_nm))) + np.square(np.mean(np.sin(psi_nm)))
        
        result = dict()
        result['gamma2_nm'] = gamma2_nm