我们从Python开源项目中,提取了以下50个代码示例,用于说明如何使用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 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
def run(self, norm_pts = None): self.exp.run_sweeps() data = {} var = {} for buff in self.exp.buffers: if self.exp.writer_to_qubit[buff.name][0] in self.qubit_names: dataset, descriptor = buff.get_data(), buff.get_descriptor() qubit_name = self.exp.writer_to_qubit[buff.name][0] if norm_pts: buff_data = normalize_data(dataset, zero_id = norm_pts[qubit_name][0], one_id = norm_pts[qubit_name][1]) else: buff_data = dataset['Data'] data[qubit_name] = self.quad_fun(buff_data) if 'Variance' in dataset.dtype.names: if self.quad in ['real', 'imag']: var[qubit_name] = self.quad_fun(dataset['Variance'])/descriptor.metadata["num_averages"] else: raise Exception('Variance of {} not available. Choose real or imag'.format(self.quad)) else: var[qubit_name] = None # Return data and variance of the mean if len(data) == 1: # if single qubit, get rid of dictionary data = list(data.values())[0] var = list(var.values())[0] return data, var
def update_references(self, frequency): # store decimated reference for mix down # phase_drift = 2j*np.pi*0.5e-6 * (abs(frequency) - 100e6) ref = np.exp(2j*np.pi * -frequency * self.time_pts[::self.d1] + 1j*self._phase, dtype=np.complex64) self.reference = ref self.reference_r = np.real(ref) self.reference_i = np.imag(ref)