Python numpy 模块，iscomplex() 实例源码

def spectrum(x,
             n0=0,
             T_s=1,
             oversample=1,
             only_positive=True):
    """
    Return the spectrum for the signal *x* calculated via FFT and the
    associated frequencies as a tuple. The *n0* parameter gives the
    index in *x* for time index 0 (*n0* = 0 means that `x[0]` is at
    time 0). The number of spectral samples returned is the next power
    of 2 greater than the length of *x* multiplied by *oversample*. If
    *only_positive*, return the spectrum only for positive frequencies
    (and raise an exception if *x* is not real).
    """
    assert oversample >= 1 and isinstance(oversample, int)
    N = nextpow2(len(x)) * 2**(oversample - 1)
    X = NP.fft.fft(x, n=N) * T_s
    f = NP.fft.fftfreq(N, d=T_s)
    if n0 != 0:
        X *= NP.exp(-1j * 2 * math.pi * NP.arange(N) * n0 / N)
    X = NP.fft.fftshift(X)
    f = NP.fft.fftshift(f)
    if only_positive:
        if any(NP.iscomplex(x)):
            raise ValueError('x is complex and only returning information for positive frequencies --- this is likely not what you want to do')
        X = X[f >= 0]
        f = f[f >= 0]
    return X, f

def plot_data(data, scroll_axis=2):
    """ Plot an image associated data.
    Currently support on 1D, 2D or 3D data.

    Parameters
    ----------
    data: array
        the data to be displayed.
    scroll_axis: int (optional, default 2)
        the scroll axis for 3d data.
    """
    # Check input parameters
    if data.ndim not in range(1, 4):
        raise ValueError("Unsupported data dimension.")

    # Deal with complex data
    if numpy.iscomplex(data).any():
        data = numpy.abs(data)

    # Create application
    app = pyqtgraph.mkQApp()

    # Create the widget
    if data.ndim == 3:
        indices = [i for i in range(3) if i != scroll_axis]
        indices = [scroll_axis] + indices
        widget = pyqtgraph.image(numpy.transpose(data, indices))
    elif data.ndim == 2:
        widget = pyqtgraph.image(data)
    else:
        widget = pyqtgraph.plot(data)

    # Run application
    app.exec_()

def average_structure(X):
    """
    Calculate an average structure from an ensemble of structures
    (i.e. X is a rank-3 tensor: X[i] is a (N,3) configuration matrix).

    @param X: m x n x 3 input vector
    @type X: numpy array

    @return: average structure
    @rtype: (n,3) numpy.array
    """
    from numpy.linalg import eigh

    B = csb.numeric.gower_matrix(X)
    v, U = eigh(B)
    if numpy.iscomplex(v).any():
        v = v.real
    if numpy.iscomplex(U).any():
        U = U.real

    indices = numpy.argsort(v)[-3:]
    v = numpy.take(v, indices, 0)
    U = numpy.take(U, indices, 1)

    x = U * numpy.sqrt(v)
    i = 0
    while is_mirror_image(x, X[0]) and i < 2:
        x[:, i] *= -1
        i += 1
    return x

def rel_entropy_true(p, q):
    """KL divergence (relative entropy) D(p||q) in bits

    Returns a scalar entropy when the input distributions p and q are
    vectors of probability masses, or returns in a row vector the
    columnwise relative entropies of the input probability matrices p and
    q"""

    if type(p) == list or type(q) == tuple:
        p = np.array(p)
    if type(q) == list or type(q) == tuple:
        q = np.array(q)

    if not p.shape == q.shape:
        raise Exception('p and q must be equal sizes',
                        'p: ' + str(p.shape),
                        'q: ' + str(q.shape))

    if (np.iscomplex(p).any() or not
        np.isfinite(p).any() or
        (p < 0).any() or
        (p > 1).any()):
        raise Exception('The probability elements of p must be real numbers'
                        'between 0 and 1.')

    if (np.iscomplex(q).any() or not
        np.isfinite(q).any() or
        (q < 0).any() or
        (q > 1).any()):
        raise Exception('The probability elements of q must be real numbers'
                        'between 0 and 1.')

    eps = math.sqrt(np.spacing(1))
    if (np.abs(np.sum(p, axis=0) - 1) > eps).any():
        raise Exception('Sum of the probability elements of p must equal 1.')
    if (np.abs(np.sum(q, axis=0) - 1) > eps).any():
        raise Exception('Sum of the probability elements of q must equal 1.')

    return sum(np.log2((p**p) / (q**p)))

def pnorm(self, p):
        r"""
        Calculates the p-norm of a vector.

        Parameters
        ----------
        p : int
            Used in computing the p-norm. This should only be set
            when calculating the pnorm of a vector is desired.

        Returns
        -------
        pn : float
            The p-norm of the vector.

        Notes
        -----
        The p-norm, which is considered a class of vector norms is defined as:

        .. math::

            ||x||_p = \sqrt[p]{|x_1|^p + |x_2|^p + \cdots + |x_n|^p} \qquad p \geq 1

        References
        ----------
        Golub, G., & Van Loan, C. (2013). Matrix computations (3rd ed.). Baltimore (MD): Johns Hopkins U.P.

        """
        if p != np.floor(p):
            p = np.floor(p)

        if np.iscomplex(p) or p < 1:
            raise ValueError('p must be at least 1 and real')

        pn = np.sum(np.absolute(self.x) ** p) ** (1. / p)

        return pn

def __str__(self):
        # create local dictionary to convert node-numbers to node-labels
        num2node_label = {num: name for name, num in self.node_label2num.items()}

        # build message to print
        msg = '------------------------\n'
        msg += '    SpicePy.Network:\n'
        msg += '------------------------\n'
        for ele, nodes, val in zip(self.names, self.nodes, self.values):
            # if val is a list --> ele is a transient source
            if isinstance(val, list):
                if self.source_type[ele] == 'pwl':
                    fmt = "{} {} {} {}(" + "{} " * (len(val[0]) - 1) + "{})\n"
                    msg += fmt.format(ele, num2node_label[nodes[0]], num2node_label[nodes[1]], self.source_type[ele], *val[0])
                else:
                    fmt = "{} {} {} {}(" + "{} " * (len(val) - 1) + "{})\n"
                    msg += fmt.format(ele, num2node_label[nodes[0]], num2node_label[nodes[1]], self.source_type[ele], *val)
            # if val is complex --> ele is a phasor
            elif np.iscomplex(val):
                msg += "{} {} {} {} {}\n".format(ele, num2node_label[nodes[0]], num2node_label[nodes[1]], np.abs(val), np.angle(val) * 180/np.pi)
            # if ele is C or L
            elif ele[0].upper() == 'C' or ele[0].upper() == 'L':
                # check if an i.c. is present and print it
                if ele in self.IC:
                    msg += "{} {} {} {} ic={}\n".format(ele, num2node_label[nodes[0]], num2node_label[nodes[1]], val, self.IC[ele])
                # otherwise...
                else:
                    msg += "{} {} {} {}\n".format(ele, num2node_label[nodes[0]], num2node_label[nodes[1]], val)
            # otherwise...general case -->  ele n+ n- val
            else:
                msg += "{} {} {} {}\n".format(ele, num2node_label[nodes[0]], num2node_label[nodes[1]], val)

        # add analysis
        msg += " ".join(self.analysis) + '\n'

        # if a plot command is present, add it
        if self.plot_cmd is not None:
            msg += self.plot_cmd + '\n'

        # add number of nodes (reference node is included) and number of branches
        msg += '------------------------\n'
        msg += '* number of nodes {}\n'.format(self.node_num + 1)
        msg += '* number of branches {}\n'.format(len(self.names))
        msg += '------------------------\n'

        return msg

def find_zero_crossing_quadratic(x1, y1, x2, y2, x3, y3, eps = 1.0):
        """ Find zero crossing using quadratic approximation along 1d line"""
        # compute coords along 1d line
        v = x2 - x1
        v = v / np.linalg.norm(v)
        if v[v!=0].shape[0] == 0:
            logging.error('Difference is 0. Probably a bug')

        t1 = 0
        t2 = (x2 - x1)[v!=0] / v[v!=0]
        t2 = t2[0]
        t3 = (x3 - x1)[v!=0] / v[v!=0]
        t3 = t3[0]

        # solve for quad approx
        x1_row = np.array([t1**2, t1, 1])
        x2_row = np.array([t2**2, t2, 1])
        x3_row = np.array([t3**2, t3, 1])
        X = np.array([x1_row, x2_row, x3_row])
        y_vec = np.array([y1, y2, y3])
        try:
            w = np.linalg.solve(X, y_vec)
        except np.linalg.LinAlgError:
            logging.error('Singular matrix. Probably a bug')
            return None

        # get positive roots
        possible_t = np.roots(w)
        t_zc = None
        for i in range(possible_t.shape[0]):
            if possible_t[i] >= 0 and possible_t[i] <= 10 and not np.iscomplex(possible_t[i]):
                t_zc = possible_t[i]

        # if no positive roots find min
        if np.abs(w[0]) < 1e-10:
            return None

        if t_zc is None:
            t_zc = -w[1] / (2 * w[0])

        if t_zc < -eps or t_zc > eps:
            return None

        x_zc = x1 + t_zc * v
        return x_zc

def run(X,y1,y2,dmax=None):

    D,N = X.shape
    y1_unique,J1 = np.unique(y1, return_inverse=True)
    ny1 = y1_unique.size

    y2_unique,J2 = np.unique(y2, return_inverse=True)
    ny2 = y2_unique.size

    Y1 = np.zeros((ny1,N))
    Y2 = np.zeros((ny2,N))
    Y1[J1,range(N)] = 1.
    Y2[J2,range(N)] = 1.

    XY2 = np.dot(X,Y2.T) # D x ny2
    XY2Y2X = np.dot(XY2,XY2.T) # D x D
    XX = np.dot(X,X.T) # D x D

    P = np.zeros((D,0))
    Sj = np.dot(X,Y1.T) #  D x ny1

    if dmax==None:
        dmax = D

    for d in range(dmax):
        if d>0: 
            invPP = np.linalg.pinv(np.dot(P.T,P))
            Sj -= np.dot(np.dot(np.dot(P,invPP),P.T),Sj) 

        C = np.dot(Sj,Sj.T) - XY2Y2X
        C = 0.5*(C+C.T)
        dd,E = scipy.linalg.eigh(C,eigvals=(D-1,D-1)) # ascending order

        assert np.isnan(dd).any()==False
        assert np.iscomplex(dd).any()==False
        #dd = dd[::-1] #
        #E = E[:,::-1]
        wj = E#E[:,0] # D x 1
        pj = np.dot(XX,wj) / np.dot(np.dot(wj.T,XX),wj) # D x 1
        P = np.hstack((P,pj.reshape((D,1)))) #  D x d


    #P = P/np.tile(np.sqrt((P**2).sum(axis=0,keepdims=True)),(D,1))
    #% They need not be orthogonal.
    return P

def run(X,y1,y2):
    # function [E,dd] = privacyLDA(X,y1,y2)
    # %max_W0 tr(W0'*C1*W0) - tr(W0'*C2*W0)
    D,N = X.shape
    y1_unique = np.unique(y1)
    #print y1_unique
    #[y1_unique,~,J1] = unique(y1);
    #ny1 = y1_unique.size

    y2_unique = np.unique(y2)
    #print y2_unique
    #[y2_unique,~,J2] = unique(y2);
    #ny2 = y2_unique.size

    C1 = np.zeros((D,D))
    C2 = np.zeros((D,D))
    mu = X.mean(axis=1).reshape((D,1))
    #print mu.shape

    for k in np.nditer(y1_unique):
        indk = np.where(y1==k)[0]        
        muk = X[:,indk].mean(axis=1).reshape((D,1))
        #muk -= np.kron(np.ones((1,len(indk))),mu)
        #%C1 = C1 + ny1*(muk-mu)*(muk-mu)';
        C1 = C1 + len(indk)*np.dot(muk-mu,(muk-mu).T)

    for k in np.nditer(y2_unique):
        indk = np.where(y2==k)[0]        
        muk = X[:,indk].mean(axis=1).reshape((D,1))
        #muk -= np.kron(np.ones((1,len(indk))),mu)
        #%C1 = C1 + ny1*(muk-mu)*(muk-mu)';
        C2 = C2 + len(indk)*np.dot(muk-mu,(muk-mu).T)

    C1 = C1 + 1e-8*np.trace(C1)*np.eye(D)# 
    C2 = C2 + 1e-8*np.trace(C2)*np.eye(D)#
    C1 = 0.5*(C1+C1.T)#;% + 1E-8*trace(C1)*eye(D); 
    C2 = 0.5*(C2+C2.T)#;% + 1E-8*trace(C2)*eye(D);


    dd,E = scipy.linalg.eigh(C1,C2) # ascending order

    #print dd.shape
    #print E.shape

    assert np.isnan(dd).any()==False
    assert np.iscomplex(dd).any()==False
    #[dd,ind] = sort(diag(dd),'descend'); 
    #print dd
    dd = dd[::-1] #
    E = E[:,::-1]

    E = E/np.tile(np.sqrt((E**2).sum(axis=0,keepdims=True)),(D,1))
    #% They need not be orthogonal.

    #print dd.shape
    #print E.shape

    return (E,dd)

def cwplot(fb_est,rx,t,fs:int,fn) -> None:
#%% time
    fg,axs = subplots(1,2,figsize=(12,6))
    ax = axs[0]
    ax.plot(t, rx.T.real)
    ax.set_xlabel('time [sec]')
    ax.set_ylabel('amplitude')
    ax.set_title('Noisy, jammed receive signal')
    #%% periodogram
    if DTPG >= (t[-1]-t[0]):
        dt = (t[-1]-t[0])/4
    else:
        dt = DTPG

    dtw = 2*dt #  seconds to window
    tstep = ceil(dt*fs)
    wind = ceil(dtw*fs);
    Nfft = zeropadfactor*wind

    f,Sraw = signal.welch(rx.ravel(), fs,
                          nperseg=wind,noverlap=tstep,nfft=Nfft,
                          return_onesided=False)

    if np.iscomplex(rx).any():
        f = np.fft.fftshift(f); Sraw = np.fft.fftshift(Sraw)

    ax=axs[1]
    ax.plot(f,Sraw,'r',label='raw signal')

    fc_est = f[Sraw.argmax()]

    #ax.set_yscale('log')
    ax.set_xlim([fc_est-200,fc_est+200])
    ax.set_xlabel('frequency [Hz]')
    ax.set_ylabel('amplitude')
    ax.legend()

    esttxt=''

    if fn is None: # simulation
        ax.axvline(ft+fb0,color='red',linestyle='--',label='true freq.')
        esttxt += f'true: {ft+fb0} Hz '

    for e in fb_est:
        ax.axvline(e,color='blue',linestyle='--',label='est. freq.')

    esttxt += ' est: ' + str(fb_est) +' Hz'

    ax.set_title(esttxt)

def cw_est(rx, fs:int, Ntone:int, method:str='esprit', usepython=False, useall=False):
    """
    estimate beat frequency using subspace frequency estimation techniques.
    This is much faster in Fortran, but to start using Python alone doesn't require compiling Fortran.

    ESPRIT and RootMUSIC are two popular subspace techniques.

    Matlab's rootmusic is a far inferior FFT-based method with very poor accuracy vs. my implementation.
    """
    assert isinstance(method,str)
    method = method.lower()

    tic = time()
    if method == 'esprit':
#%% ESPRIT
        if rx.ndim == 2:
            assert usepython,'Fortran not yet configured for multi-pulse case'
            Ntone *= 2


        if usepython or (Sc is None and Sr is None):
            print('Python ESPRIT')
            fb_est, sigma = esprit(rx, Ntone, Nblockest, fs)
        elif np.iscomplex(rx).any():
            print('Fortran complex64 ESPRIT')
            fb_est, sigma = Sc.subspace.esprit(rx,Ntone,Nblockest,fs)
        else: # real signal
            print('Fortran float32 ESPRIT')
            fb_est, sigma = Sr.subspace.esprit(rx,Ntone,Nblockest,fs)

        fb_est = abs(fb_est)
#%% ROOTMUSIC
    elif method == 'rootmusic':
        fb_est, sigma = rootmusic(rx,Ntone,Nblockest,fs)
    else:
        raise ValueError(f'unknown estimation method: {method}')
    print(f'computed via {method} in {time()-tic:.1f} seconds.')
#%% improvised process for CW only without notch filter
    # assumes first two results have largest singular values (from SVD)
    if not useall:
        i = sigma > 0.001 # arbitrary
        fb_est = fb_est[i]
        sigma  = sigma[i]

#        if fb_est.size>1:
#            ii = np.argpartition(sigma, Ntone-1)[:Ntone-1]
#            fb_est = fb_est[ii]
#            sigma = sigma[ii]


    return fb_est, sigma