Python scipy.angle函数代码示例

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)