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

def phormants(x, Fs):
    N = len(x)
    w = numpy.hamming(N)

    # Apply window and high pass filter.
    x1 = x * w   
    x1 = lfilter([1], [1., 0.63], x1)

    # Get LPC.    
    ncoeff = 2 + Fs / 1000
    A, e, k = lpc(x1, ncoeff)    
    #A, e, k = lpc(x1, 8)

    # Get roots.
    rts = numpy.roots(A)
    rts = [r for r in rts if numpy.imag(r) >= 0]

    # Get angles.
    angz = numpy.arctan2(numpy.imag(rts), numpy.real(rts))

    # Get frequencies.    
    frqs = sorted(angz * (Fs / (2 * math.pi)))

    return frqs

def mdst(x, odd=True):
    """ Calculate modified discrete sine transform of input signal

    Parameters
    ----------
    X : array_like
        The input signal
    odd : boolean, optional
        Switch to oddly stacked transform. Defaults to :code:`True`.

    Returns
    -------
    out : array_like
        The output signal

    """
    return -1 * numpy.imag(cmdct(x, odd=odd)) * numpy.sqrt(2)

def get_phases(self):
        sizeimg = np.real(self.imgfft).shape
        mag = np.zeros(sizeimg)
        for x in range(sizeimg[0]):
            for y in range(sizeimg[1]):
                mag[x][y] = np.arctan2(np.real(self.imgfft[x][y]), np.imag(self.imgfft[x][y]))
        rpic = MyImage(mag)
        rpic.limit(1)
        return rpic    


#      int my = y-output.height/2;
#      int mx = x-output.width/2;
#      float angle = atan2(my, mx) - HALF_PI ;
#      float radius = sqrt(mx*mx+my*my) / factor;
#      float ix = map(angle,-PI,PI,input.width,0);
#      float iy = map(radius,0,height,0,input.height);
#      int inputIndex = int(ix) + int(iy) * input.width;
#      int outputIndex = x + y * output.width;
#      if (inputIndex <= input.pixels.length-1) {
#        output.pixels[outputIndex] = input.pixels[inputIndex];

def fftDf( df , part = "abs") :
   #Handle series or DataFrame
   if type(df) == pd.Series :
      df = pd.DataFrame(df)
      ise = True
   else :
      ise = False
   res = pd.DataFrame( index = np.fft.rfftfreq( df.index.size, d = dx( df ) ) )
   for col in df.columns :
      if part == "abs" :
         res["FFT_"+col] = np.abs( np.fft.rfft(df[col]) )  / (0.5*df.index.size)
      elif part == "real" :
         res["FFT_"+col] = np.real( np.fft.rfft(df[col]) ) / (0.5*df.index.size)
      elif part == "imag" :
         res["FFT_"+col] = np.imag( np.fft.rfft(df[col]) ) / (0.5*df.index.size)
   if ise :
      return res.iloc[:,0]
   else :
      return res

def test_psi(adjcube):
    """Tests retrieval of the wave functions and eigenvalues.
    """
    from pydft.bases.fourier import psi, O, H
    cell = adjcube
    V = QHO(cell)
    W = W4(cell)
    Ns = W.shape[1]
    Psi, epsilon = psi(V, W, cell, forceR=False)

    #Make sure that the eigenvalues are real.
    assert np.sum(np.imag(epsilon)) < 1e-13

    checkI = np.dot(Psi.conjugate().T, O(Psi, cell))
    assert abs(np.sum(np.diag(checkI))-Ns) < 1e-13 # Should be the identity
    assert np.abs(np.sum(checkI)-Ns) < 1e-13

    checkD = np.dot(Psi.conjugate().T, H(V, Psi, cell))
    diagsum = np.sum(np.diag(checkD))
    assert np.abs(np.sum(checkD)-diagsum) < 1e-12 # Should be diagonal

    # Should match the diagonal elements of previous matrix
    assert np.allclose(np.diag(checkD), epsilon)

def fcn_ComputeFrequencyResponse(self,f,sig,mur,a,x0,y0,z0,X,Y,Z):
        """Compute Single Frequency Response at (X,Y,Z)"""

        m = self.m
        orient = self.orient
        xtx = self.xtx
        ytx = self.ytx
        ztx = self.ztx

        chi = fcn_ComputeExcitation_FEM(f,sig,mur,a)
        Hpx,Hpy,Hpz = fcn_ComputePrimary(m,orient,xtx,ytx,ztx,x0,y0,z0)

        mx = 4*np.pi*a**3*chi*Hpx/3
        my = 4*np.pi*a**3*chi*Hpy/3
        mz = 4*np.pi*a**3*chi*Hpz/3
        R = np.sqrt((X-x0)**2 + (Y-y0)**2 + (Z-z0)**2)

        Hx = (1/(4*np.pi))*(3*(X-x0)*(mx*(X-x0) + my*(Y-y0) + mz*(Z-z0))/R**5 - mx/R**3)
        Hy = (1/(4*np.pi))*(3*(Y-y0)*(mx*(X-x0) + my*(Y-y0) + mz*(Z-z0))/R**5 - my/R**3)
        Hz = (1/(4*np.pi))*(3*(Z-z0)*(mx*(X-x0) + my*(Y-y0) + mz*(Z-z0))/R**5 - mz/R**3)
        Habs = np.sqrt(np.real(Hx)**2 + np.real(Hy)**2 + np.real(Hz)**2) + 1j*np.sqrt(np.imag(Hx)**2 + np.imag(Hy)**2 + np.imag(Hz)**2)

        return Hx, Hy, Hz, Habs

def genSpectra(time,dipole,signal):

    fw, frequency = pade(time,dipole)
    fw_sig, frequency = pade(time,signal,alternate=True)

    fw_re = np.real(fw)
    fw_im = np.imag(fw)
    fw_abs = fw_re**2 + fw_im**2

    #spectra = (fw_re*17.32)/(np.pi*field*damp_const)
    #spectra = (fw_re*17.32*514.220652)/(np.pi*field*damp_const)
    #numerator = np.imag((fw*np.conjugate(fw_sig)))
    numerator = np.imag(fw_abs*np.conjugate(fw_sig))
    #numerator = np.abs((fw*np.conjugate(fw_sig)))
    #numerator = np.abs(fw)
    denominator = np.real(np.conjugate(fw_sig)*fw_sig)
    #denominator = 1.0 
    spectra = ((4.0*27.21138602*2*frequency*np.pi*(numerator))/(3.0*137.036*denominator))
    spectra *= 1.0/100.0
    #plt.plot(frequency*27.2114,fourier)
    #plt.show()
    return frequency, spectra

def histogram_plot(data, sfreq, toffset, bins, log_scale, title):
    """Plot a histogram of the data for a given bin size."""
    print("histogram")

    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.hist(numpy.real(data), bins,
            log=log_scale, histtype='bar', color=['green'])
    ax.hold(True)
    ax.hist(numpy.imag(data), bins,
            log=log_scale, histtype='bar', color=['blue'])
    ax.grid(True)
    ax.set_xlabel('adc value')
    ax.set_ylabel('frequency')
    ax.set_title(title)
    ax.hold(False)

    return fig

def __init__(self,jet,kernels,k,x,y,pt,subpixel):
        self.jet = jet
        self.kernels = kernels
        self.k = k
        self.x = x
        self.y = y

        re = np.real(jet)
        im = np.imag(jet)

        self.mag = np.sqrt(re*re + im*im)
        self.phase = np.arctan2(re,im)

        if subpixel:
            d = np.array([[pt.X()-x],[pt.Y()-y]])
            comp = np.dot(self.k,d)
            self.phase -= comp.flatten()
            self.jet = self.mag*np.exp(1.0j*self.phase)

def __init__(self, qubit_names, quad="real"):
        super(PulseCalibration, self).__init__()
        self.qubit_names = qubit_names if isinstance(qubit_names, list) else [qubit_names]
        self.qubit     = [QubitFactory(qubit_name) for qubit_name in qubit_names] if isinstance(qubit_names, list) else QubitFactory(qubit_names)
        self.filename   = 'None'
        self.exp        = None
        self.axis_descriptor = None
        self.cw_mode    = False
        self.saved_settings = config.load_meas_file(config.meas_file)
        self.settings = deepcopy(self.saved_settings) #make a copy for used during calibration
        self.quad = quad
        if quad == "real":
            self.quad_fun = np.real
        elif quad == "imag":
            self.quad_fun = np.imag
        elif quad == "amp":
            self.quad_fun = np.abs
        elif quad == "phase":
            self.quad_fun = np.angle
        else:
            raise ValueError('Quadrature to calibrate must be one of ("real", "imag", "amp", "phase").')
        self.plot       = self.init_plot()

def fit_photon_number(xdata, ydata, params):
    ''' Fit number of measurement photons before a Ramsey. See McClure et al., Phys. Rev. App. 2016
    input params:
    1 - cavity decay rate kappa (MHz)
    2 - detuning Delta (MHz)
    3 - dispersive shift 2Chi (MHz)
    4 - Ramsey decay time T2* (us)
    5 - exp(-t_meas/T1) (us), only if starting from |1> (to include relaxation during the 1st msm't)
    6 - initial qubit state (0/1)
    '''
    params = [2*np.pi*p for p in params[:3]] + params[3:] # convert to angular frequencies
    def model_0(t, pa, pb):
        return (-np.imag(np.exp(-(1/params[3]+params[1]*1j)*t + (pa-pb*params[2]*(1-np.exp(-((params[0] + params[2]*1j)*t)))/(params[0]+params[2]*1j))*1j)))
    def model(t, pa, pb):
        return  params[4]*model_0(t, pa, pb) + (1-params[4])*model_0(t, pa+np.pi, pb) if params[5] == 1  else model_0(t, pa, pb)
    popt, pcov = curve_fit(model, xdata, ydata, p0 = [0, 1])
    perr = np.sqrt(np.diag(pcov))
    finer_delays = np.linspace(np.min(xdata), np.max(xdata), 4*len(xdata))
    fit_curve = model(finer_delays, *popt)
    return popt[1], perr[1], (finer_delays, fit_curve)

def make_layout(self):
        self.lay = QtWidgets.QHBoxLayout()
        self.lay.setContentsMargins(0, 0, 0, 0)
        self.real = FloatSpinBox(label=self.labeltext,
                                 min=self.minimum,
                                 max=self.maximum,
                                 increment=self.singleStep,
                                 log_increment=self.log_increment,
                                 halflife_seconds=self.halflife_seconds,
                                 decimals=self.decimals)
        self.imag = FloatSpinBox(label=self.labeltext,
                                 min=self.minimum,
                                 max=self.maximum,
                                 increment=self.singleStep,
                                 log_increment=self.log_increment,
                                 halflife_seconds=self.halflife_seconds,
                                 decimals=self.decimals)
        self.real.value_changed.connect(self.value_changed)
        self.lay.addWidget(self.real)
        self.label = QtWidgets.QLabel(" + j")
        self.lay.addWidget(self.label)
        self.imag.value_changed.connect(self.value_changed)
        self.lay.addWidget(self.imag)
        self.setLayout(self.lay)
        self.setFocusPolicy(QtCore.Qt.ClickFocus)

def set_value(self, obj, value):
        """
        the master's setter writes its value to the slave lists
        """
        real, complex = [], []
        for v in value:
            # separate real from complex values
            if np.imag(v) == 0:
                real.append(v.real)
            else:
                complex.append(v)
        # avoid calling setup twice
        with obj.do_setup:
            setattr(obj, 'complex_' + self.name, complex)
            setattr(obj, 'real_' + self.name, real)
        # this property should have call_setup=True, such that obj._setup()
        # is called automatically after this function

def plot_waveforms(waveforms, figTitle=''):
    channels = waveforms.keys()
    # plot
    plots = []
    for (ct, chan) in enumerate(channels):
        fig = bk.figure(title=figTitle + repr(chan),
                        plot_width=800,
                        plot_height=350,
                        y_range=[-1.05, 1.05],
                        x_axis_label=u'Time (?s)')
        fig.background_fill_color = config.plotBackground
        if config.gridColor:
            fig.xgrid.grid_line_color = config.gridColor
            fig.ygrid.grid_line_color = config.gridColor
        waveformToPlot = waveforms[chan]
        xpts = np.linspace(0, len(waveformToPlot) / chan.phys_chan.sampling_rate
                           / 1e-6, len(waveformToPlot))
        fig.line(xpts, np.real(waveformToPlot), color='red')
        fig.line(xpts, np.imag(waveformToPlot), color='blue')
        plots.append(fig)
    bk.show(column(*plots))

def merge_waveform(n, chAB, chAm1, chAm2, chBm1, chBm2):
    '''
    Builds packed I and Q waveforms from the nth mini LL, merging in marker data.
    '''
    wfAB = np.array([], dtype=np.complex)
    for entry in chAB['linkList'][n % len(chAB['linkList'])]:
        if not entry.isTimeAmp:
            wfAB = np.append(wfAB, chAB['wfLib'][entry.key])
        else:
            wfAB = np.append(wfAB, chAB['wfLib'][entry.key][0] *
                             np.ones(entry.length * entry.repeat))

    wfAm1 = marker_waveform(chAm1['linkList'][n % len(chAm1['linkList'])],
                            chAm1['wfLib'])
    wfAm2 = marker_waveform(chAm2['linkList'][n % len(chAm2['linkList'])],
                            chAm2['wfLib'])
    wfBm1 = marker_waveform(chBm1['linkList'][n % len(chBm1['linkList'])],
                            chBm1['wfLib'])
    wfBm2 = marker_waveform(chBm2['linkList'][n % len(chBm2['linkList'])],
                            chBm2['wfLib'])

    wfA = pack_waveform(np.real(wfAB), wfAm1, wfAm2)
    wfB = pack_waveform(np.imag(wfAB), wfBm1, wfBm2)

    return wfA, wfB

def check(value, value_list, difference):
    n = True
    if len(value_list) == 0:
        value_list.append(value)
    else:
        for x in value_list:
            if np.abs(np.real(x) - np.real(value)) < difference and \
               np.abs(np.imag(x) - np.imag(value)) < difference:
                   n = False
            else:
                pass
        if n == True:
            value_list.append(value)
    return value_list

# This function converts a list of lists into a numpy array. It only takes the
# list of lists as input, and returns the array as output. If the lists inside 
# the list are of unequal lengths, it fills up the lines with None so that all 
# lines in the output array are of equal length.
# Example input:
# a = [[1,3,4], [2,1], [2,3,4,7]]
# Output:
# array([[1, 3, 4, None],
#        [2, 1, None, None],
#        [2, 3, 4, 7]], dtype=object)

def csvd(arr):
    """
    Do the complex SVD of a 2D array, returning real valued U, S, VT

    http://stemblab.github.io/complex-svd/
    """
    C_r = arr.real
    C_i = arr.imag
    block_x = C_r.shape[0]
    block_y = C_r.shape[1]
    K = np.zeros((2 * block_x, 2 * block_y))
    # Upper left
    K[:block_x, :block_y] = C_r
    # Lower left
    K[:block_x, block_y:] = C_i
    # Upper right
    K[block_x:, :block_y] = -C_i
    # Lower right
    K[block_x:, block_y:] = C_r
    return svd(K, full_matrices=False)

def csvd(arr):
    """
    Do the complex SVD of a 2D array, returning real valued U, S, VT

    http://stemblab.github.io/complex-svd/
    """
    C_r = arr.real
    C_i = arr.imag
    block_x = C_r.shape[0]
    block_y = C_r.shape[1]
    K = np.zeros((2 * block_x, 2 * block_y))
    # Upper left
    K[:block_x, :block_y] = C_r
    # Lower left
    K[:block_x, block_y:] = C_i
    # Upper right
    K[block_x:, :block_y] = -C_i
    # Lower right
    K[block_x:, block_y:] = C_r
    return svd(K, full_matrices=False)

def fft_test2(self):
        axis = str(self.axis_combobox.currentText())

        if axis.startswith('a'):
            normal_para = 16384.0
        elif axis.startswith('g'):
            normal_para = 131.0
        signal =( self.raw_data[axis] - self.bias_dict[axis])/ normal_para

        n = signal.size # Number of data points
        dx = 0.007 # Sampling period (in meters)
        Fk = np.fft.fft(signal) # Fourier coefficients (divided by n)
        nu = np.fft.fftfreq(n,dx) # Natural frequencies
        #Fk = np.fft.fftshift(Fk) # Shift zero freq to center
        #nu = np.fft.fftshift(nu) # Shift zero freq to center
        f, ax = plt.subplots(3,1,sharex=True)
        ax[0].plot(nu, np.real(Fk)) # Plot Cosine terms
        ax[0].set_ylabel(r'$Re[F_k]$', size = 'x-large')
        ax[1].plot(nu, np.imag(Fk)) # Plot Sine terms
        ax[1].set_ylabel(r'$Im[F_k]$', size = 'x-large')
        ax[2].plot(nu, np.absolute(Fk)**2) # Plot spectral power
        ax[2].set_ylabel(r'$\vert F_k \vert ^2$', size = 'x-large')
        ax[2].set_xlabel(r'$\widetilde{\nu}$', size = 'x-large')
        plt.title(axis)
        plt.show()

def estimate_pair(self, ts1, ts2):
        """

        Returns
        -------
        ts : array-like, shape(1, n_samples)
            Estimated iPLV time series.

        avg : float
            Average iPLV.


        Notes
        -----
        Called from :mod:`dyfunconn.tvfcgs.tvfcg`.
        """
        n_samples = len(ts1)

        ts_plv = np.exp(1j * (ts1 - ts2))
        avg_plv = np.abs(np.imag(np.sum((ts_plv))) / float(n_samples))

        return np.imag(ts_plv), avg_plv

def edge_phase():
    """calculate edge phase"""
    se = plane.UniformPlane(L=8, W=8, js=(0, 8 * 7), E=0, t=1, U=0, phase=.2 * 2 * np.pi, parasite=.1)

    E1, psi1l, psi1r = eigenbasis(se, 1)
    idx = np.argsort(np.real(E1))
    E1 = E1[idx]
    psi1l = psi1l[:, idx]
    psi1r = psi1r[:, idx]

    res = np.zeros((64, ))
    idxs = se.edge_indices(dw=1, dl=1)
    print(idxs)
    s = len(idxs)

    for i in range(s):
        res += np.array([np.arctan2(np.real(psi1r[idxs[i], j] / psi1r[idxs[(i + 1) % s], j]), np.imag(psi1r[idxs[i], j] / psi1l[idxs[(i + 1) % s], j])) for j in np.arange(64)])

    plt.plot(np.real(E1), res / (2 * np.pi), '-o')
    Emin = np.min(np.real(E1))
    Emax = np.max(np.real(E1))
    for i in range(-10, 1, 1):
        plt.plot([Emin, Emax], [i, i])
        plt.plot([Emin, Emax], [-i, -i])
    plt.show()

def BB(Y,index_PQ, index_P, n_PQ, n_P):
    case_number, _ = np.shape(Y)
    Y_p = Y.copy()
    B_p = np.zeros((n_P,n_P))
    B_pp = np.zeros((n_PQ,n_PQ))
    #--------------------------------------------------
    for i in xrange(case_number):
        Y_p[i][i] = complex(0,0)
        for j in xrange(case_number):
            if i != j:
                Y_p[i][i] -= Y_p[i][j]

    B = np.imag(Y_p)
    for i in xrange(n_P):
        for j in xrange(0, n_P):
            B_p[i][j] = B[index_P[i]][index_P[j]]
    #--------------------------------------------------
    for i in xrange(0, n_PQ):
        for j in xrange(0, n_PQ):
            B_pp[i][j] = B[index_PQ[i]][index_PQ[j]]

    return B_p, B_pp

# A.M Van Amerongen-----------------------------------------------------------------------------------------------------

def vals2coeffs2(vals):
    """Map function values at Chebyshev points of 2nd kind to
    first-kind Chebyshev polynomial coefficients"""
    n = vals.size
    if n <= 1:
        coeffs = vals
        return coeffs
    tmp = np.append( vals[::-1], vals[1:-1] )
    if np.isreal(vals).all():
        coeffs = ifft(tmp)
        coeffs = np.real(coeffs)
    elif np.isreal( 1j*vals ).all():
        coeffs = ifft(np.imag(tmp))
        coeffs = 1j * np.real(coeffs)
    else:
        coeffs = ifft(tmp)
    coeffs = coeffs[:n]
    coeffs[1:n-1] = 2*coeffs[1:n-1]
    return coeffs

def coeffs2vals2(coeffs):
    """Map first-kind Chebyshev polynomial coefficients to
    function values at Chebyshev points of 2nd kind"""
    n = coeffs.size
    if n <= 1:
        vals = coeffs
        return vals
    coeffs = coeffs.copy()
    coeffs[1:n-1] = .5 * coeffs[1:n-1]
    tmp = np.append( coeffs, coeffs[n-2:0:-1] )
    if np.isreal(coeffs).all():
        vals = fft(tmp)
        vals = np.real(vals)
    elif np.isreal(1j*coeffs).all():
        vals = fft(np.imag(tmp))
        vals = 1j * np.real(vals)
    else:
        vals = fft(tmp)
    vals = vals[n-1::-1]
    return vals

def _plot_samples(self, signal, ax, mag, real, imag, rms, noise=True):
        if mag:
            ax.plot(signal.mag, label='Mag')
        if real:
            ax.plot(np.real(signal), label='Real')
        if imag:
            ax.plot(np.imag(signal), label='Imag')
        if rms:
            ax.axhline(signal.rms, label='RMS', linestyle='--')
        if noise:
            noise_est = self.result.carrier_info.noise / np.sqrt(len(signal))
            ax.axhline(noise_est, label='Noise', linestyle='--', color='g')
        ax.legend()
        ax.set_xlabel('Sample')
        ax.set_ylabel('Value')
        # ax2 = ax.twiny()
        # ax2.set_xlim(0, len(signal) / self.sample_rate * 1e3)
        # ax2.set_xlabel('Time (ms)')
        ax.grid()

def interpolate_slice(f3d, rot, pfac=2, size=None):
    nhalf = f3d.shape[0] / 2
    if size is None:
        phalf = nhalf
    else:
        phalf = size / 2
    qot = rot * pfac  # Scaling!
    px, py, pz = np.meshgrid(np.arange(-phalf, phalf), np.arange(-phalf, phalf), 0)
    pr = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
    pcoords = np.vstack([px.reshape(-1), py.reshape(-1), pz.reshape(-1)])
    mcoords = qot.T.dot(pcoords)
    mcoords = mcoords[:, pr.reshape(-1) < nhalf]
    pvals = map_coordinates(np.real(f3d), mcoords, order=1, mode="wrap") + \
             1j * map_coordinates(np.imag(f3d), mcoords, order=1, mode="wrap")
    pslice = np.zeros(pr.shape, dtype=np.complex)
    pslice[pr < nhalf] = pvals
    return pslice

def complex_quadrature(func, a, b, **kwargs):
    """
    wraps the scipy qaudpack routines to handle complex valued functions

    :param func: callable
    :param a: lower limit
    :param b: upper limit
    :param kwargs: kwargs for func
    :return:
    """

    def real_func(x):
        return np.real(func(x))

    def imag_func(x):
        return np.imag(func(x))

    real_integral = integrate.quad(real_func, a, b, **kwargs)
    imag_integral = integrate.quad(imag_func, a, b, **kwargs)

    return real_integral[0] + 1j * imag_integral[0], real_integral[1] + imag_integral[1]

def plot(self):
        """ Plot a realisation of the signal waveform """
        Y=self.rvs()
        Y_processed=linear_transform(Y,self.preprocessing_method)
        N,L=Y_processed.shape

        if ((L==3) or (L==1)):
            n_vect=np.arange(N)/self.Fe
            for l in range(L):
                plt.plot(n_vect,Y_processed[:,l],label="signal %d" %l)
            plt.xlabel("Time")
            plt.ylabel("Signal")
            plt.legend()

        if L==2:
            z=Y_processed[:,0]+1j*Y_processed[:,1]
            plt.plot(np.real(z),np.imag(z))
            plt.xlabel("Real Part")
            plt.ylabel("Imag Part")

def make_node(self, a, s=None):
        a = T.as_tensor_variable(a)
        if a.ndim < 3:
            raise TypeError('%s: input must have dimension >= 3,  with ' %
                            self.__class__.__name__ +
                            'first dimension batches and last real/imag parts')

        if s is None:
            s = a.shape[1:-1]
            s = T.set_subtensor(s[-1], (s[-1] - 1) * 2)
            s = T.as_tensor_variable(s)
        else:
            s = T.as_tensor_variable(s)
            if (not s.dtype.startswith('int')) and \
               (not s.dtype.startswith('uint')):
                raise TypeError('%s: length of the transformed axis must be'
                                ' of type integer' % self.__class__.__name__)
        return gof.Apply(self, [a, s], [self.output_type(a)()])

def add_scal_vec(self, val, vec):
        """
        Perform in-place addition of a vector times a scalar.

        Parameters
        ----------
        val : int or float
            scalar.
        vec : <Vector>
            this vector times val is added to self.
        """
        if self._vector_info._under_complex_step:
            r_val = np.real(val)
            i_val = np.imag(val)
            for set_name, data in iteritems(self._data):
                data += r_val * vec._data[set_name] + i_val * vec._imag_data[set_name]
            for set_name, data in iteritems(self._imag_data):
                data += i_val * vec._data[set_name] + r_val * vec._imag_data[set_name]
        else:
            for set_name, data in iteritems(self._data):
                data += val * vec._data[set_name]

def get_batch(batch_size):
    samples = np.zeros([batch_size, sample_length])
    frequencies = [set()] * batch_size
    ffts = np.zeros([batch_size, fft_size])

    for i in range(batch_size):
        num_sources = np.random.randint(min_sources, max_sources + 1)
        for source_idx in range(num_sources):
            frequency, sample = generate_sample()
            samples[i] += sample
            frequencies[i].add(frequency)

        samples[i] /= float(num_sources)

        fft = np.fft.rfft(samples[i], norm="ortho")
        fft = np.real(fft)**2 + np.imag(fft)**2

        fft *= fft_norm

        ffts[i] = fft

    return frequencies, samples, ffts

def get_imag_part(self):
        r = MyImage(np.imag(self.imgfft))  
        r.limit(1)
        return r
    # Correlate functions

def get_magnitude(self):
        sizeimg = np.real(self.imgfft).shape
        mag = np.zeros(sizeimg)
        for x in range(sizeimg[0]):
            for y in range(sizeimg[1]):
                mag[x][y] = np.sqrt(np.real(self.imgfft[x][y])**2 + np.imag(self.imgfft[x][y])**2)
        rpic = MyImage(mag)
        rpic.limit(1)
        return rpic

def draw_fft(self):
        if len(self.points) < 1:
            return
        pts = map(lambda p: p[1] - self.offset, self.points)

        out = numpy.fft.rfft(pts)
        c = len(out)

        norm = 0
        for i in range(c/2):
            norm += numpy.real(out[i])**2 + numpy.imag(out[i])**2

        norm = math.sqrt(norm)
        if norm <= 0:
            return

        for i in range(1, SignalKPlot.NUM_X_DIV):
            x = float(i) / SignalKPlot.NUM_X_DIV
            glRasterPos2d(x, .95)
            period = 3/math.exp(x) # incorrect!!
            SignalKPlot.drawputs(str(period))

        glPushMatrix()
        glBegin(GL_LINE_STRIP)
        for i in range(c/2):
            glVertex2d(float(i) * 2 / (c-2), abs(out[i]) / norm)
        glEnd()
        glPopMatrix()

def slidingFFT( se, T  ,  n = 1 , tStart = None , preSample = False , nHarmo = 5 , kind = abs , phase = None) :
   """
   Harmonic analysis on a sliding windows
   se : Series to analyse
   T : Period
   tStart : start _xAxis
   n : size of the sliding windows in period.
   reSample : reSample the signal so that a period correspond to a integer number of time step
   nHarmo : number of harmonics to return
   kind : module, real,  imaginary part, as a function (abs, np.imag, np.real ...)
   phase : phase shift (for instance to extract in-phase with cos or sin)
   """

   if (type(se) == pd.DataFrame) :
      if len(se.columns) == 1 : se = se.iloc[:,0]

   nWin = int(0.5 + n*T / dx(se) )
   #ReSample to get round number of time step per period
   if preSample :
      new = reSample( se, dt = n*T / (nWin) )
   else :
      new = se
   signal = new.values[:]
   #Allocate results
   res = np.zeros( (new.shape[0] , nHarmo ) )
   for iWin in range(new.shape[0] - nWin) :
      sig = signal[ iWin : iWin+nWin  ]  #windows
      fft = np.fft.fft( sig )            #FTT
      if phase !=None  :                 #Phase shift
         fft *= np.exp( 1j* ( 2*pi*(iWin*1./nWin) + phase ))
      fftp = kind( fft )       #Take module, real or imaginary part
      spectre = 2*fftp/(nWin)  #Scale
      for ih in range(nHarmo):
         res[iWin, ih] = spectre[ih*n]
         if ih == 0 : res[iWin, ih] /= 2.0
         #if ih == 0 : res[iWin, ih] = 2.0
   return pd.DataFrame( data = res , index = new.index , columns = map( lambda x : "Harmo {:} ({:})".format(x , se.name  ) , range(nHarmo)  ) )

def test_E_real(adjcube):
    """Tests that the result of the calculation is real.
    """
    from pydft.bases.fourier import E
    from numpy.matlib import randn
    cell = adjcube

    #Single columns of random complex data
    W = np.array(randn(np.prod(cell.S), 4) + 1j*randn(np.prod(cell.S), 4))
    #Setup a harmonic oscillator potential
    V = QHO(cell)
    En = E(V, W, cell, forceR=False)

    assert np.imag(En) < 1e-14

def test_IJ(adjcube):
    """Tests the I and J operators."""
    from pydft.bases.fourier import I, J
    #This also tests accessing the geometry via the global variable.
    Sprod = np.prod(adjcube.S)
    for i in range(10):
        v = np.random.random(size=Sprod)
        #Our v is real; but due to round-off problems, there will be
        #tiny imaginary values. Chop them off.
        it = J(I(v))
        if abs(np.max(np.imag(it))) < 1e-14:
            it = np.real(it)
        assert np.allclose(it, v)

def test_LLinv(adjcube):
    """Tests L and its inverse.
    """
    from pydft.bases.fourier import L, Linv
    Sprod = np.prod(adjcube.S)
    for i in range(10):
        v = np.random.random(size=Sprod)
        #Our v is real; but due to round-off problems, there will be
        #tiny imaginary values. Chop them off. We only keep the last
        #N-1 components because the 0 component is NaN.
        it = Linv(L(v))[1:]
        if abs(np.max(np.imag(it))) < 1e-14:
            it = np.real(it)
        assert np.allclose(it, v[1:])

def plotResponseFEM(Ax,fi,f,H,Comp):

    FS = 20

    xTicks = (np.logspace(np.log(np.min(f)),np.log(np.max(f)),9))
    Ylim = np.array([np.min(np.real(H)),np.max(np.real(H))])

    Ax.grid('both', linestyle='-', linewidth=0.8, color=[0.8, 0.8, 0.8])
    Ax.semilogx(f,0*f,color='k',linewidth=2)
    Ax.semilogx(f,np.real(H),color='k',linewidth=4,label="Real")
    Ax.semilogx(f,np.imag(H),color='k',linewidth=4,ls='--',label="Imaginary")
    Ax.semilogx(np.array([fi,fi]),1.1*Ylim,linewidth=3,color='r')
    Ax.set_xbound(np.min(f),np.max(f))
    Ax.set_ybound(1.1*Ylim)
    Ax.set_xlabel('Frequency [Hz]',fontsize=FS+2)
    Ax.tick_params(labelsize=FS-2)
    Ax.yaxis.set_major_formatter(FormatStrFormatter('%.1e'))

    if Comp == 'x':
        Ax.set_ylabel('$\mathbf{Hx}$ [A/m]',fontsize=FS+4,labelpad=-5)
        Ax.set_title('$\mathbf{Hx}$ Response at $\mathbf{Rx}$',fontsize=FS+6)
    elif Comp == 'y':
        Ax.set_ylabel('$\mathbf{Hy}$ [A/m]',fontsize=FS+4,labelpad=-5)
        Ax.set_title('$\mathbf{Hy}$ Response at $\mathbf{Rx}$',fontsize=FS+6)
    elif Comp == 'z':
        Ax.set_ylabel('$\mathbf{Hz}$ [A/m]',fontsize=FS+4,labelpad=-5)
        Ax.set_title('$\mathbf{Hz}$ Response at $\mathbf{Rx}$',fontsize=FS+6)
    elif Comp == 'abs':
        Ax.set_ylabel('$\mathbf{|H|}$ [A/m]',fontsize=FS+4,labelpad=-5)
        Ax.set_title('$\mathbf{|H|}$ Response at $\mathbf{Rx}$',fontsize=FS+6)


    if np.max(np.real(H[-1])) > 0.:
        handles, labels = Ax.get_legend_handles_labels()
        Ax.legend(handles, labels, loc='upper left', fontsize=FS)
    elif np.max(np.real(H[-1])) < 0.:
        handles, labels = Ax.get_legend_handles_labels()
        Ax.legend(handles, labels, loc='lower left', fontsize=FS)

    return Ax

def plot_InducedCurrent_FD(self,Ax,Is,fi):

        FS = 20

        R = self.R
        L = self.L

        Imax = np.max(-np.real(Is))

        f = np.logspace(0,8,101)


        Ax.grid('both', linestyle='-', linewidth=0.8, color=[0.8, 0.8, 0.8])
        Ax.semilogx(f,-np.real(Is),color='k',linewidth=4,label="$I_{Re}$")
        Ax.semilogx(f,-np.imag(Is),color='k',ls='--',linewidth=4,label="$I_{Im}$")
        Ax.semilogx(fi*np.array([1.,1.]),np.array([0,1.1*Imax]),color='r',ls='-',linewidth=3)
        handles, labels = Ax.get_legend_handles_labels()
        Ax.legend(handles, labels, loc='upper left', fontsize=FS)

        Ax.set_xlabel('Frequency [Hz]',fontsize=FS+2)
        Ax.set_ylabel('$\mathbf{- \, I_s (\omega)}$ [A]',fontsize=FS+2,labelpad=-10)
        Ax.set_title('Frequency Response',fontsize=FS)
        Ax.set_ybound(0,1.1*Imax)
        Ax.tick_params(labelsize=FS-2)
        Ax.yaxis.set_major_formatter(FormatStrFormatter('%.1e'))

        #R_str    = '{:.3e}'.format(R)
        #L_str    = '{:.3e}'.format(L)
        #f_str    = '{:.3e}'.format(fi)
        #EMF_str  = '{:.2e}j'.format(EMFi.imag)
        #I_str    = '{:.2e} - {:.2e}j'.format(float(np.real(Isi)),np.abs(float(np.imag(Isi))))

        #Ax.text(1.4,1.01*Imax,'$R$ = '+R_str+' $\Omega$',fontsize=FS)
        #Ax.text(1.4,0.94*Imax,'$L$ = '+L_str+' H',fontsize=FS)
        #Ax.text(1.4,0.87*Imax,'$f$ = '+f_str+' Hz',fontsize=FS,color='r')
        #Ax.text(1.4,0.8*Imax,'$V$ = '+EMF_str+' V',fontsize=FS,color='r')
        #Ax.text(1.4,0.73*Imax,'$I_s$ = '+I_str+' A',fontsize=FS,color='r')

        return Ax

def test_real(self):
        y = np.random.rand(10,)
        assert_array_equal(0, np.imag(y))

def test_cmplx(self):
        y = np.random.rand(10,)+1j*np.random.rand(10,)
        assert_array_equal(y.imag, np.imag(y))

def test_complex_bad2(self):
        with np.errstate(divide='ignore', invalid='ignore'):
            v = 1 + 1j
            v += np.array(-1+1.j)/0.
        vals = nan_to_num(v)
        assert_all(np.isfinite(vals))
        # Fixme
        #assert_all(vals.imag > 1e10)  and assert_all(np.isfinite(vals))
        # !! This is actually (unexpectedly) positive
        # !! inf.  Comment out for now, and see if it
        # !! changes
        #assert_all(vals.real < -1e10) and assert_all(np.isfinite(vals))

def voltage_plot(data, sfreq, toffset, log_scale, title):
    """Plot the real and imaginary voltage from IQ data."""

    print("voltage")

    t_axis = numpy.arange(0, len(data)) / sfreq + toffset

    fig = plt.figure()
    ax0 = fig.add_subplot(2, 1, 1)
    ax0.plot(t_axis, data.real)
    ax0.grid(True)
    maxr = numpy.max(data.real)
    minr = numpy.min(data.real)

    if minr == 0.0 and maxr == 0.0:
        minr = -1.0
        maxr = 1.0

    ax0.axis([t_axis[0], t_axis[len(t_axis) - 1], minr, maxr])
    ax0.set_ylabel('I sample value (A/D units)')

    ax1 = fig.add_subplot(2, 1, 2)
    ax1.plot(t_axis, data.imag)
    ax1.grid(True)
    maxi = numpy.max(data.imag)
    mini = numpy.min(data.imag)

    if mini == 0.0 and maxi == 0.0:
        mini = -1.0
        maxi = 1.0

    ax1.axis([t_axis[0], t_axis[len(t_axis) - 1], mini, maxi])

    ax1.set_xlabel('time (seconds)')

    ax1.set_ylabel('Q sample value (A/D units)')
    ax1.set_title(title)

    return fig

def iq_plot(data, toffset, log_scale, title):
    """Plot an IQ circle from the data in linear or log scale."""
    print("iq")

    if log_scale:
        rx_raster_r = numpy.sign(
            data.real) * numpy.log10(numpy.abs(data.real) + 1E-30) / numpy.log10(2.)
        rx_raster_i = numpy.sign(
            data.imag) * numpy.log10(numpy.abs(data.imag) + 1E-30) / numpy.log10(2.)
    else:
        data *= 1.0 / 32768.0
        rx_raster_r = data.real
        rx_raster_i = data.imag

    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(rx_raster_r, rx_raster_i, '.')

    axmx = numpy.max([numpy.max(rx_raster_r), numpy.max(rx_raster_i)])

    ax.axis([-axmx, axmx, -axmx, axmx])
    ax.grid(True)
    ax.set_xlabel('I')
    ax.set_ylabel('Q')
    ax.set_title(title)

    return fig

def show(self,*args,**kwargs):
        print self.data.shape
        tiles = []
        w,h,n = self.data.shape
        for i in range(n):
            mat = self.data[:,:,i]
            tiles.append(pv.Image(mat.real))
            tiles.append(pv.Image(mat.imag))
        mont = pv.ImageMontage(tiles,layout=(8,10),tileSize=(w,h))
        mont.show(*args,**kwargs)

def test_gabor1(self):
        ilog = None # pv.ImageLog(name="GaborTest1")

        bank = FilterBank(tile_size=(128,128))
        kernels = createGaborKernels()
        for wavelet in kernels:
            bank.addFilter(wavelet)

        for i in range(len(bank.filters)):
            filter = np.fft.ifft2(bank.filters[i])
            if ilog:
                ilog.log(pv.Image(np.fft.fftshift(filter.real)),label="Filter_RE_%d"%i)
                ilog.log(pv.Image(np.fft.fftshift(filter.imag)),label="Filter_IM_%d"%i)


        for im in self.test_images[:1]:
            if ilog:
                ilog.log(im,label="ORIG")
            results = bank.convolve(im)
            #print "RShape",results.shape[2]
            if ilog:
                for i in range(results.shape[2]):
                    ilog.log(pv.Image(results[:,:,i].real),label="CONV_RE_%d"%i)
                    ilog.log(pv.Image(results[:,:,i].imag),label="CONV_IM_%d"%i)
        if ilog:
            ilog.show()

def generate_fake_data(alpha, phi, sigma, N = 5000, plot=False):

    N_samples = 256
    data_start = 3
    data_length = 100
    gnd_mean = np.array([alpha*np.cos(phi), alpha*np.sin(phi)])
    ex_mean = np.array([alpha*np.cos(phi + np.pi), alpha*np.sin(phi + np.pi)])
    gndIQ = np.vectorize(complex)(np.random.normal(gnd_mean[0], sigma, N),
                                 np.random.normal(gnd_mean[1], sigma, N))
    exIQ = np.vectorize(complex)(np.random.normal(ex_mean[0], sigma, N),
                                 np.random.normal(ex_mean[1], sigma, N))
    gnd = np.zeros((N_samples, N), dtype=np.complex128)
    ex = np.zeros((N_samples, N), dtype=np.complex128)
    for idx, x in enumerate(zip(gndIQ, exIQ)):
        gnd[data_start:data_start+data_length, idx] = x[0]
        ex[data_start:data_start+data_length, idx] = x[1]

    gnd += sigma/50 * (np.random.randn(N_samples, N) + 1j * np.random.randn(N_samples, N))
    ex += sigma/50 * (np.random.randn(N_samples, N) + 1j * np.random.randn(N_samples, N))

    if plot:
        plt.figure()
        plt.plot(np.real(gndIQ), np.imag(gndIQ), 'b.')
        plt.plot(np.real(exIQ), np.imag(exIQ), 'r.')
        plt.draw()
        plt.show()

        plt.figure()
        plt.plot(np.real(gnd[:,15]), 'b.')
        plt.plot(np.real(ex[:,15]), 'r.')
        plt.draw()
        plt.show()
    return gnd, ex