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