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
python类imag()的实例源码
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]
