The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
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 Example 2
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) Example 3
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]; Example 4
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 Example 5
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) Example 6
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 Example 7
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 Example 8
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 Example 9
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) Example 10
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() Example 11
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) Example 12
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) Example 13
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 Example 14
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)) Example 15
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 Example 16
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) Example 17
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) Example 18
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) Example 19
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() Example 20
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 Example 21
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() Example 22
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----------------------------------------------------------------------------------------------------- Example 23
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 Example 24
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 Example 25
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() Example 26
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 Example 27
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] Example 28
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") Example 29
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)()]) Example 30
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] Example 31
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 Example 32
def get_imag_part(self):
r = MyImage(np.imag(self.imgfft))
r.limit(1)
return r
# Correlate functions Example 33
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 Example 34
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() Example 35
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) ) ) Example 36
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 Example 37
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) Example 38
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:]) Example 39
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 Example 40
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 Example 41
def test_real(self):
y = np.random.rand(10,)
assert_array_equal(0, np.imag(y)) Example 42
def test_cmplx(self):
y = np.random.rand(10,)+1j*np.random.rand(10,)
assert_array_equal(y.imag, np.imag(y)) Example 43
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)) Example 44
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 Example 45
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 Example 46
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) Example 47
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() Example 48
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 Example 49
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 Example 50
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) 