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 generate_mtf(pixel_aperture=1, azimuth=0, num_samples=128):
''' generates the 1D diffraction-limited MTF for a given pixel size and azimuth.
Args:
pixel_aperture (`float`): aperture of the pixel, in microns. Pixel is
assumed to be square.
azimuth (`float`): azimuth to retrieve the MTF at, in degrees.
num_samples (`int`): number of samples in the output array.
Returns:
`tuple` containing:
`numpy.ndarray`: array of units, in cy/mm.
`numpy.ndarray`: array of MTF values (rel. 1.0).
Notes:
Azimuth is not actually implemented yet.
'''
pitch_unit = pixel_aperture / 1e3
normalized_frequencies = np.linspace(0, 2, num_samples)
otf = np.sinc(normalized_frequencies)
mtf = np.abs(otf)
return normalized_frequencies / pitch_unit, mtf Example 2
def smPhiDP(phiDP, ran):
# smooth phiDP field and take derivative
# calculate lanczos filter weights
numRan = ran.shape[0]
numK = 31
fc = 0.015
kt = np.linspace(-(numK-1)/2, (numK-1)/2, numK)
w = np.sinc(2.*kt*fc)*(2.*fc)*np.sinc(kt/(numK/2))
#smoothPhiDP = convolve1d(phiDP, w, axis=1, mode='constant', cval=-999.)
smoothPhiDP = conv(phiDP, w)
#smoothPhiDP = np.ma.masked_where(smoothPhiDP==-999., smoothPhiDP)
return smoothPhiDP
# function for estimating kdp
#---------------------------------- Example 3
def n_even_fcn(f, o, w, l):
"""Even case."""
# Variables :
k = np.array(range(0, int(l) + 1, 1)) + 0.5
b = np.zeros(k.shape)
# # Run Loop :
for s in range(0, len(f), 2):
m = (o[s + 1] - o[s]) / (f[s + 1] - f[s])
b1 = o[s] - m * f[s]
b = b + (m / (4 * np.pi * np.pi) * (np.cos(2 * np.pi * k * f[
s + 1]) - np.cos(2 * np.pi * k * f[s])) / (
k * k)) * abs(np.square(w[round((s + 1) / 2)]))
b = b + (f[s + 1] * (m * f[s + 1] + b1) * np.sinc(2 * k * f[
s + 1]) - f[s] * (m * f[s] + b1) * np.sinc(2 * k * f[s])) * abs(
np.square(w[round((s + 1) / 2)]))
a = (np.square(w[0])) * 4 * b
h = 0.5 * np.concatenate((np.flipud(a), a))
return h Example 4
def NevenFcn(F, M, W, L): # N is even
# Variables :
k = np.array(range(0, int(L) + 1, 1)) + 0.5
b = np.zeros(k.shape)
# # Run Loop :
for s in range(0, len(F), 2):
m = (M[s + 1] - M[s]) / (F[s + 1] - F[s])
b1 = M[s] - m * F[s]
b = b + (m / (4 * np.pi * np.pi) * (np.cos(2 * np.pi * k * F[
s + 1]) - np.cos(2 * np.pi * k * F[s])) / (
k * k)) * abs(np.square(W[round((s + 1) / 2)]))
b = b + (F[s + 1] * (m * F[s + 1] + b1) * np.sinc(2 * k * F[
s + 1]) - F[s] * (m * F[s] + b1) * np.sinc(2 * k * F[s])) * abs(
np.square(W[round((s + 1) / 2)]))
a = (np.square(W[0])) * 4 * b
h = 0.5 * np.concatenate((np.flipud(a), a))
return h
####################################################################
# - Filt the signal :
#################################################################### Example 5
def lpfls(N,wp,ws,W):
M = (N-1)/2
nq = np.arange(0,2*M+1)
nb = np.arange(0,M+1)
q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
b[0] = wp/np.pi
q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
b = b.transpose()
Q1 = ln.toeplitz(q[0:M+1])
Q2 = ln.hankel(q[0:M+1],q[M:])
Q = Q1+Q2
a = ln.solve(Q,b)
h = list(nq)
for i in nb:
h[i] = 0.5*a[M-i]
h[N-1-i] = h[i]
h[M] = 2*h[M]
hmax = max(np.absolute(h))
for i in nq:
h[i] = (8191/hmax)*h[i]
return h Example 6
def bpfls(N,ws1,wp1,wp2,ws2,W):
M = (N-1)/2
nq = np.arange(0,2*M+1)
nb = np.arange(0,M+1)
q = W*np.sinc(nq) - (W*ws2/np.pi) * np.sinc(nq* (ws2/np.pi)) + (wp2/np.pi) * np.sinc(nq*(wp2/np.pi)) - (wp1/np.pi) * np.sinc(nq*(wp1/np.pi)) + (W*ws1/np.pi) * np.sinc(nq*(ws1/np.pi))
b = (wp2/np.pi)*np.sinc((wp2/np.pi)*nb) - (wp1/np.pi)*np.sinc((wp1/np.pi)*nb)
b[0] = wp2/np.pi - wp1/np.pi
q[0] = W - W*ws2/np.pi + wp2/np.pi - wp1/np.pi + W*ws1/np.pi # since sin(pi*n)/pi*n = 1, not 0
b = b.transpose()
Q1 = ln.toeplitz(q[0:M+1])
Q2 = ln.hankel(q[0:M+1],q[M:])
Q = Q1+Q2
a = ln.solve(Q,b)
h = list(nq)
for i in nb:
h[i] = 0.5*a[M-i]
h[N-1-i] = h[i]
h[M] = 2*h[M]
hmax = max(np.absolute(h))
for i in nq:
h[i] = (8191/hmax)*h[i]
return h Example 7
def hpfls(N,ws,wp,W):
M = (N-1)/2
nq = np.arange(0,2*M+1)
nb = np.arange(0,M+1)
b = 1 - (wp/np.pi)* np.sinc(nb * wp/np.pi)
b[0] = 1- wp/np.pi
q = 1 - (wp/np.pi)* np.sinc(nq * wp/np.pi) + W * (ws/np.pi) * np.sinc(nq * ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
q[0] = b[0] + W* ws/np.pi
b = b.transpose()
Q1 = ln.toeplitz(q[0:M+1])
Q2 = ln.hankel(q[0:M+1],q[M:])
Q = Q1+Q2
a = ln.solve(Q,b)
h = list(nq)
for i in nb:
h[i] = 0.5*a[M-i]
h[N-1-i] = h[i]
h[M] = 2*h[M]
hmax = max(np.absolute(h))
for i in nq:
h[i] = (8191/hmax)*h[i]
return h Example 8
def lanczosSubPixKernel( subPixShift, kernelShape=3, lobes=None ):
"""
Generate a kernel suitable for ni.convolve to subpixally shift an image.
"""
kernelShape = np.array( [kernelShape], dtype='int' )
if kernelShape.ndim == 1: # make it 2-D
kernelShape = np.array( [kernelShape[0], kernelShape[0]], dtype='int' )
if lobes is None:
lobes = (kernelShape[0]+1)/2
x_range = np.arange(-kernelShape[1]/2,kernelShape[1]/2)+1.0-subPixShift[1]
x_range = ( 2.0 / kernelShape[1] ) * x_range
y_range = np.arange(-kernelShape[1]/2,kernelShape[0]/2)+1.0-subPixShift[0]
y_range = ( 2.0 /kernelShape[0] ) * y_range
[xmesh,ymesh] = np.meshgrid( x_range, y_range )
lanczos_filt = np.sinc(xmesh * lobes) * np.sinc(xmesh) * np.sinc(ymesh * lobes) * np.sinc(ymesh)
lanczos_filt = lanczos_filt / np.sum(lanczos_filt) # Normalize filter output
return lanczos_filt Example 9
def sincinterp(x):
"""
Sinc interpolation for computation of fractional transformations.
As appears in :
-https://github.com/audiolabs/frft/
----------
Args:
f : (array) Complex valued input array
a : (float) Alpha factor
Returns:
ret : (array) Real valued synthesised data
"""
N = len(x)
y = np.zeros(2 * N - 1, dtype=x.dtype)
y[:2 * N:2] = x
xint = fftconvolve( y[:2 * N], np.sinc(np.arange(-(2 * N - 3), (2 * N - 2)).T / 2),)
return xint[2 * N - 3: -2 * N + 3] Example 10
def collect_pixel(self, pixel_i, frame_i, k,j,i):
#print pixel_i, k,j,i
t0 = time.time()
#px_data = np.random.rand()
#px_data = t0 - self.prev_px
x_hw = self.stage.settings.x_position.read_from_hardware(send_signal=False)
y_hw = self.stage.settings.y_position.read_from_hardware(send_signal=False)
theta = np.pi/10. * frame_i
x = x_hw*np.cos(theta) - y_hw*np.sin(theta)
y = x_hw*np.sin(theta) + y_hw*np.cos(theta)
px_data = np.sinc((x-50)*0.05)**2 * np.sinc(0.05*(y-50))**2 #+ 0.05*np.random.random()
#px_data = (x-xhw)**2 + ( y-yhw)**2
#if px_data > 1:
# print('hw', x, xhw, y, yhw)
self.display_image_map[k,j,i] = px_data
if self.settings['save_h5']:
self.test_data[frame_i, k,j,i] = px_data
time.sleep(self.settings['pixel_time'])
#self.prev_px = t0 Example 11
def design_windowed_sinc_lpf(fc, bw):
N = Filter.get_filter_length_from_bandwidth(bw)
# Compute sinc filter impulse response
h = np.sinc(2 * fc * (np.arange(N) - (N - 1) / 2.))
# We use blackman window function
w = np.blackman(N)
# Multiply sinc filter with window function
h = h * w
# Normalize to get unity gain
h_unity = h / np.sum(h)
return h_unity Example 12
def intensitiesFFIntraAtom(nq, dq, partNrs, nameList, ffDict):
""" uses atomistic form factors """
dIntensity = np.zeros((nq, 2), dtype=float)
partInt = np.zeros((nq, len(partNrs) + 1))
qList = np.zeros(nq)
qList[:] = [float(i * dq) for i in range(nq)]
dIntensity[:, 0] = qList[:]
partInt[:, 0] = qList[:]
# for j in range(0,len(dIntegrand)):
# r=dIntegrand[j,0]
# sinc=j0(q*r)
k = 0
formFacProd = np.zeros((nq, len(partNrs) + 1))
# partNrsProd=getPartNrsProd(partNrs)
# print "partNrsProd ", partNrsProd
for i in range(nq):
for k in range(1, len(partNrs) + 1):
# print k
formFacProd[i, k] = ff.fiveGaussian(
ffDict[nameList[k - 1]], qList[i]) ** 2
partInt[i, k] += partNrs[k - 1] * formFacProd[i, k]
for i in range(nq):
dIntensity[i, 1] = partInt[i, 1:].sum()
return partInt, dIntensity Example 13
def iterate_l1(L, alpha, arg, beta, K, N, rr):
oversample_ratio = (1.0 * K / N)
for l1 in range(-L, L + 1):
alf = alpha[abs(l1)] * 1.0
if l1 < 0:
alf = numpy.conj(alf)
# r1 = numpy.sinc(1.0*(arg+1.0*l1*beta)/(1.0*K/N))
input_array = (arg + 1.0 * l1 * beta) / oversample_ratio
r1 = dirichlet(input_array.astype(numpy.float32))
rr = iterate_sum(rr, alf, r1)
return rr Example 14
def nufft_r(om, N, J, K, alpha, beta):
'''
equation (30) of Fessler's paper
'''
def iterate_sum(rr, alf, r1):
rr = rr + alf * r1
return rr
def iterate_l1(L, alpha, arg, beta, K, N, rr):
oversample_ratio = (1.0 * K / N)
import time
t0=time.time()
for l1 in range(-L, L + 1):
alf = alpha[abs(l1)] * 1.0
# if l1 < 0:
# alf = numpy.conj(alf)
# r1 = numpy.sinc(1.0*(arg+1.0*l1*beta)/(1.0*K/N))
input_array = (arg + 1.0 * l1 * beta) / oversample_ratio
r1 = dirichlet(input_array)
rr = iterate_sum(rr, alf, r1)
return rr
M = numpy.size(om) # 1D size
gam = 2.0 * numpy.pi / (K * 1.0)
nufft_offset0 = nufft_offset(om, J, K) # om/gam - nufft_offset , [M,1]
dk = 1.0 * om / gam - nufft_offset0 # om/gam - nufft_offset , [M,1]
arg = outer_sum(-numpy.arange(1, J + 1) * 1.0, dk)
L = numpy.size(alpha) - 1
# print('alpha',alpha)
rr = numpy.zeros((J, M), dtype=numpy.float32)
rr = iterate_l1(L, alpha, arg, beta, K, N, rr)
return (rr, arg) Example 15
def analytic_ft(self, unit_x, unit_y):
''' Analytic fourier transform of a pixel aperture.
Args:
unit_x (`numpy.ndarray`): sample points in x axis.
unit_y (`numpy.ndarray`): sample points in y axis.
Returns:
`numpy.ndarray`: 2D numpy array containing the analytic fourier transform.
'''
xq, yq = np.meshgrid(unit_x, unit_y)
return (sinc(xq * self.size_x / 1e3) *
sinc(yq * self.size_y / 1e3)).astype(config.precision) Example 16
def make_bin_weights(self):
erb_max = hz2erb(self.sr/2.0)
ngrid = 1000
erb_grid = np.arange(ngrid) * erb_max / (ngrid - 1)
hz_grid = (np.exp(erb_grid / 9.26) - 1) / 0.00437
resp = np.zeros((ngrid, self.n_bins))
for b in range(self.n_bins):
w = self.widths[b]
r = (2.0 * w + 1.0) / self.sr * (hz_grid - self.hz_freqs[b])
resp[:,b] = np.square(np.sinc(r)+ 0.5 * np.sinc(r + 1.0) + 0.5 * np.sinc(r - 1.0));
self.weights = np.dot(linalg.pinv(resp), np.ones((ngrid,1))) Example 17
def autocorrelation_function(self, r):
"""compute the real space autocorrelation function for the Gaussian random field model
"""
# compute the cut-level parameter beta
beta = np.sqrt(2) * erfinv(2 * (1-self.frac_volume) - 1)
# the covariance of the GRF
acf_psi = (np.exp(-r/self.corr_length) * (1 + r / self.corr_length)
* np.sinc(2*r/self.repeat_distance))
# integral discretization. henning says: the resolution 1e-2 is ad hoc, test required,
# the integrand has a (integrable) singularity for t=1 and acf_psi = 1, so an adaptive
# discretization seems preferable -> TODO
dt = 1e-2
t = np.arange(0, 1, dt)
# the gridded integrand, via change of integration variable
# compared to the wp-2 docu, to enable array-based computation
t_gridded, acf_psi_gridded = np.meshgrid(t, acf_psi)
integrand_gridded = (acf_psi_gridded / np.sqrt(1 - (t_gridded * acf_psi_gridded)**2)
* np.exp( - beta**2 / (1 + t_gridded * acf_psi_gridded)))
acf = 1.0 / (2 * np.pi) * np.trapz(integrand_gridded, x=t_gridded)
return acf
# ft not known analytically: deligate
# def ft_autocorrelation_function(self, k): Example 18
def autocorrelation_function(self, r):
"""compute the real space autocorrelation function for the Teubner Strey model
"""
acf = self.corr_func_at_origin * np.exp(-r/self.corr_length) * np.sinc(2*r/self.repeat_distance)
return acf Example 19
def n_odd_fcn(f, o, w, l):
"""Odd case."""
# Variables :
b0 = 0
m = np.array(range(int(l + 1)))
k = m[1:len(m)]
b = np.zeros(k.shape)
# Run Loop :
for s in range(0, len(f), 2):
m = (o[s + 1] - o[s]) / (f[s + 1] - f[s])
b1 = o[s] - m * f[s]
b0 = b0 + (b1 * (f[s + 1] - f[s]) + m / 2 * (
f[s + 1] * f[s + 1] - f[s] * f[s])) * abs(
np.square(w[round((s + 1) / 2)]))
b = b + (m / (4 * np.pi * np.pi) * (
np.cos(2 * np.pi * k * f[s + 1]) - np.cos(2 * np.pi * k * f[s])
) / (k * k)) * abs(np.square(w[round((s + 1) / 2)]))
b = b + (f[s + 1] * (m * f[s + 1] + b1) * np.sinc(2 * k * f[
s + 1]) - f[s] * (m * f[s] + b1) * np.sinc(2 * k * f[s])) * abs(
np.square(w[round((s + 1) / 2)]))
b = np.insert(b, 0, b0)
a = (np.square(w[0])) * 4 * b
a[0] = a[0] / 2
aud = np.flipud(a[1:len(a)]) / 2
a2 = np.insert(aud, len(aud), a[0])
h = np.concatenate((a2, a[1:] / 2))
return h Example 20
def NoddFcn(F, M, W, L): # N is odd
# Variables :
b0 = 0
m = np.array(range(int(L + 1)))
k = m[1:len(m)]
b = np.zeros(k.shape)
# Run Loop :
for s in range(0, len(F), 2):
m = (M[s + 1] - M[s]) / (F[s + 1] - F[s])
b1 = M[s] - m * F[s]
b0 = b0 + (b1 * (F[s + 1] - F[s]) + m / 2 * (
F[s + 1] * F[s + 1] - F[s] * F[s])) * abs(
np.square(W[round((s + 1) / 2)]))
b = b + (m / (4 * np.pi * np.pi) * (
np.cos(2 * np.pi * k * F[s + 1]) - np.cos(2 * np.pi * k * F[s])
) / (k * k)) * abs(np.square(W[round((s + 1) / 2)]))
b = b + (F[s + 1] * (m * F[s + 1] + b1) * np.sinc(2 * k * F[
s + 1]) - F[s] * (m * F[s] + b1) * np.sinc(2 * k * F[s])) * abs(
np.square(W[round((s + 1) / 2)]))
b = np.insert(b, 0, b0)
a = (np.square(W[0])) * 4 * b
a[0] = a[0] / 2
aud = np.flipud(a[1:len(a)]) / 2
a2 = np.insert(aud, len(aud), a[0])
h = np.concatenate((a2, a[1:] / 2))
return h
# Even case Example 21
def lpfls2notch(N,wp,ws,wn1,wn2,W):
M = (N-1)/2
nq = np.arange(0,2*M+1)
nb = np.arange(0,M+1)
q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
b = np.asmatrix(b)
b = b.transpose()
Q1 = ln.toeplitz(q[0:M+1])
Q2 = ln.hankel(q[0:M+1],q[M:])
Q = Q1+Q2
G1 = np.cos(wn1*nb)
G2 = np.cos(wn2*nb)
G = np.matrix([G1,G2])
d = np.array([0,0])
d = np.asmatrix(d)
d = d.transpose()
c = np.asmatrix(ln.solve(Q,b))
mu = ln.solve(G*ln.inv(Q)*G.transpose(),G*c - d)
a = c - ln.solve(Q,G.transpose()*mu)
h = np.zeros(N)
for i in nb:
h[i] = 0.5*a[M-i]
h[N-1-i] = h[i]
h[M] = 2*h[M]
hmax = max(np.absolute(h))
for i in nq:
h[i] = (8191/hmax)*h[i]
return h Example 22
def lpfls1notch(N,wp,ws,wn1,W):
M = (N-1)/2
nq = np.arange(0,2*M+1)
nb = np.arange(0,M+1)
q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
b = np.asmatrix(b)
b = b.transpose()
Q1 = ln.toeplitz(q[0:M+1])
Q2 = ln.hankel(q[0:M+1],q[M:])
Q = Q1+Q2
G1 = np.cos(wn1*nb)
G = np.matrix([G1])
d = np.array([0])
d = np.asmatrix(d)
c = np.asmatrix(ln.solve(Q,b))
mu = ln.solve(G*ln.inv(Q)*G.transpose(),G*c - d)
a = c - ln.solve(Q,G.transpose()*mu)
h = np.zeros(N)
for i in nb:
h[i] = 0.5*a[M-i]
h[N-1-i] = h[i]
h[M] = 2*h[M]
hmax = max(np.absolute(h))
for i in nq:
h[i] = (8191/hmax)*h[i]
return h Example 23
def rcosfir(beta, sps, span=None):
"""Generates a raised cosine FIR filter.
:param beta: shape of the raised cosine filter (0-1)
:param sps: number of samples per symbol
:param span: length of the filter in symbols (None => automatic selection)
>>> import arlpy
>>> rc = arlpy.comms.rcosfir(0.25, 6)
>>> bb = arlpy.comms.modulate(arlpy.comms.random_data(100), arlpy.comms.psk())
>>> pb = arlpy.comms.upconvert(bb, 6, 27000, 18000, rc)
"""
if beta < 0 or beta > 1:
raise ValueError('Beta must be between 0 and 1')
if span is None:
# from http://www.commsys.isy.liu.se/TSKS04/lectures/3/MichaelZoltowski_SquareRootRaisedCosine.pdf
# since this recommendation is for root raised cosine filter, it is conservative for a raised cosine filter
span = 33-int(44*beta) if beta < 0.68 else 4
delay = int(span*sps/2)
t = _np.arange(-delay, delay+1, dtype=_np.float)/sps
denom = 1 - (2*beta*t)**2
eps = _np.finfo(float).eps
idx1 = _np.nonzero(_np.abs(denom) > _sqrt(eps))
b = _np.full_like(t, beta*_sin(_pi/(2*beta))/(2*sps))
b[idx1] = _np.sinc(t[idx1]) * _cos(_pi*beta*t[idx1])/denom[idx1] / sps
b /= _sqrt(_np.sum(b**2))
return b Example 24
def lanczos(x, freq, radius=3):
l = 0.5*freq * np.sinc(0.5*freq*x) * np.sinc(x/radius)
l[(x < -radius) | (x > radius)] = 0.0
return l Example 25
def pollData(self):
"""Poll for new data. This method sleeps in order to ensure
that self.pollSize observations are generated at a realistic rate.
"""
# figure time between polls using exponentially weighted moving average
curTime = time.time()
if self.pollDelay >= 0.0:
self.pollDelay = 0.8*self.pollDelay + 0.2*(curTime - self.lastPollTime)
else:
self.pollDelay = 0.0
self.lastPollTime = curTime - self.shift
sleepTime = np.max((0.0, self.shift - self.pollDelay))
self.lastPollTime = curTime + sleepTime
time.sleep(sleepTime)
with self.lock:
# generate some random data
data = np.random.uniform(-self.scale.value, self.scale.value,
size=(self.pollSize,self.nChan))
erpEnd = self.erpStart.value +\
self.erpSpeed.value *\
data.shape[0]/float(self.sampRate)
erp = np.linspace(self.erpStart.value, erpEnd, data.shape[0])
erp = np.repeat(erp, data.shape[1]).reshape((-1,data.shape[1]))
erp = erp * 0.5*(np.arange(data.shape[1])+1.0)
erp = np.sinc(erp)
data += erp
self.erpStart.value = erpEnd
return data Example 26
def lanczos(n, radius=3):
taps = np.linspace(-radius,radius, n)
return np.sinc(taps/radius) Example 27
def sinc(n, radius=3, freq=1.0):
taps = np.linspace(-radius,radius, n)
return freq * np.sinc(freq*taps) Example 28
def initImpulseResponse(self, window):
if self.bandType == 'allpass':
self.impulseResponse = windows.kroneckerDelta(self.order+1)
elif self.bandType == 'allstop':
self.impulseResponse = np.zeros_like(window)
elif self.bandType == 'lowpass':
hightaps = self.high*self.taps
self.impulseResponse = self.high*np.sinc(hightaps) * window
elif self.bandType == 'highpass':
lowtaps = self.low*self.taps
self.impulseResponse = (-self.low*np.sinc(lowtaps) * window +
windows.kroneckerDelta(self.order+1))
elif self.bandType == 'bandpass':
lowtaps = self.low*self.taps
hightaps = self.high*self.taps
self.impulseResponse = (self.high*np.sinc(hightaps) -
self.low*np.sinc(lowtaps)) * window
elif self.bandType == 'bandstop':
lowtaps = self.low*self.taps
hightaps = self.high*self.taps
self.impulseResponse = ((self.high*np.sinc(hightaps) -
self.low*np.sinc(lowtaps)) * window +
windows.kroneckerDelta(self.order+1))
else:
raise Exception('Invalid bandType: ' + str(self.bandType))
self.impulseResponse = self.impulseResponse.astype(self.dtype, copy=False) Example 29
def BandpassFilter(lowFreq, highFreq, sampRate=1.0, order=None, filtType='butter', **kwargs):
filtType = filtType.lower()
if filtType in ('butter', 'cheby1', 'cheby2', 'ellip', 'bessel'):
if order is None: order = 3
return BandpassFilterIIR(lowFreq=lowFreq, highFreq=highFreq,
sampRate=sampRate, order=order, filtType=filtType, **kwargs)
elif filtType in ('lanczos', 'sinc-blackman', 'sinc-hamming', 'sinc-hann'):
if order is None: order = 20
return BandpassFilterFIR(lowFreq=lowFreq, highFreq=highFreq,
sampRate=sampRate, order=order, filtType=filtType, **kwargs)
else:
raise Exception('Invalid filter type: ' + str(filtType) + '.') Example 30
def rootRaisedCosine(t):
bit_period = 1/BIT_FREQUENCY
if (t== bit_period/(4*BETA)):
return (BETA/(np.pi*np.sqrt(2*bit_period)) * \
((np.pi + 2)*np.sin(np.pi/(4*BETA)) + (np.pi - 2)*np.cos(np.pi/(4*BETA))))
else:
return (4 * BETA / np.pi / np.sqrt(bit_period) * \
(np.cos((1 + BETA) * np.pi * t / bit_period) + \
(1 - BETA) * np.pi / (4 * BETA) * np.sinc((1-BETA)*t/bit_period)) / \
(1 - (4*BETA*t/bit_period)**2)) Example 31
def rootRaisedCosine(t):
bit_period = 1/BIT_FREQUENCY
if (t== bit_period/(4*BETA)):
return (BETA/(np.pi*np.sqrt(2*bit_period)) * \
((np.pi + 2)*np.sin(np.pi/(4*BETA)) + (np.pi - 2)*np.cos(np.pi/(4*BETA))))
else:
return (4 * BETA / np.pi / np.sqrt(bit_period) * \
(np.cos((1 + BETA) * np.pi * t / bit_period) + \
(1 - BETA) * np.pi / (4 * BETA) * np.sinc((1-BETA)*t/bit_period)) / \
(1 - (4*BETA*t/bit_period)**2)) Example 32
def get_CML(self, q, t):
"""
Calculate C, M, L forming the elements of T matrix
:param q: a shifted coordinate grid
:param t: time
:return: tuple C, M, L
"""
assert q is self.x_plus or q is self.x_minus, \
"the shifted coordinate (q) must be either x_plus or x_minus"
# get the difference of adiabatic potential curves
Vg_minus_Ve = (self._Vg_plus_Ve_x_plus if q is self.x_plus else self._Vg_minus_Ve_x_minus)
Veg = self.Veg(q, t)
D = Veg**2 + 0.25*Vg_minus_Ve**2
np.sqrt(D, out=D)
S = np.sinc(D * self.dt / np.pi)
S *= self.dt
C = D * self.dt
np.cos(C, out=C)
M = S * Vg_minus_Ve
M *= 0.5
L = S * Veg
return C, M, L Example 33
def ef_cascade(screenpos, i, nletters):
v = np.array([0, -1])
d = lambda t : 1 if t < 0 else abs(np.sinc(t) / (1 + t**4))
return lambda t: screenpos + v * 400 * d(t - 0.15 * i) Example 34
def collect_pixel(self, pixel_i, k,j,i):
#print pixel_i, k,j,i
t0 = time.time()
#px_data = np.random.rand()
#px_data = t0 - self.prev_px
x0,y0 = self.pos
x_set = self.stage.settings['x_position']
y_set = self.stage.settings['y_position']
x_hw = self.stage.settings.x_position.read_from_hardware(send_signal=False)
y_hw = self.stage.settings.y_position.read_from_hardware(send_signal=False)
if np.abs(x_hw - x0) > 1:
self.log.debug('='*60)
self.log.debug('pos {} {}'.format(x0, y0))
self.log.debug('settings {} {}'.format(x_set, y_set))
self.log.debug('hw {} {}'.format(x_hw, y_hw))
self.log.debug('settings value delta {} {}'.format(x_set-x0, y_set-y0))
self.log.debug('read_hw value delta {} {}'.format(x_hw-x0, y_hw-y0))
self.log.debug('='*60)
x = x_hw
y = y_hw
px_data = np.sinc((x-50)*0.05)**2 * np.sinc(0.05*(y-50))**2 #+ 0.05*np.random.random()
#px_data = (x-xhw)**2 + ( y-yhw)**2
#if px_data > 1:
# print('hw', x, xhw, y, yhw)
self.display_image_map[k,j,i] = px_data
if self.settings['save_h5']:
self.test_data[k,j,i] = px_data
time.sleep(self.settings['pixel_time'])
#self.prev_px = t0 Example 35
def Sinc(freq=440, amp=1.0, offset=0):
"""Makes a Sinc function.
freq: float frequency in Hz
amp: float amplitude, 1.0 is nominal max
offset: float phase offset in radians
returns: Sinusoid object
"""
return Sinusoid(freq, amp, offset, func=np.sinc) Example 36
def intensitiesFFFaster(nq, dq, dIntegrand, keys, ffDict, dr):
""" uses atomistic form factors """
dIntensity = np.zeros((nq, 2), dtype=float)
partInt = np.zeros((nq, len(dIntegrand[0])))
nameList = [k.split(",") for k in keys]
# print "nameList=",nameList
qList = np.zeros(nq)
qList[:] = [float(i * dq) for i in range(nq)]
dIntensity[:, 0] = qList[:]
partInt[:, 0] = qList[:]
rList = dIntegrand[:, 0]
# for j in range(0,len(dIntegrand)):
# r=dIntegrand[j,0]
# sinc=j0(q*r)
formFacProd = np.zeros((nq, len(dIntegrand[0])))
for i in range(nq):
sincList = np.sinc(rList * qList[i] / math.pi) * dr
for k in range(1, len(dIntegrand[0])):
# print k
formFacProd[i, k] = ff.fiveGaussian(ffDict[nameList[k - 1][0]], qList[i])\
* ff.fiveGaussian(ffDict[nameList[k - 1][1]], qList[i])
partInt[i, k] += (sincList[:] * dIntegrand[:, k]
).sum() * formFacProd[i, k]
for i in range(nq):
dIntensity[i, 1] = partInt[i, 1:].sum()
return partInt, dIntensity Example 37
def dirichlet(x):
return numpy.sinc(x) Example 38
def dirichlet(x):
return numpy.sinc(x) Example 39
def initLanczos(self, filtOrder):
self.filtOrder = filtOrder
if self.filtOrder % 2 != 0:
raise Exception('Invalid filtOrder: ' + str(self.filtOrder) +
' Must be an even integer.')
radius = self.filtOrder // 2
win = np.sinc(np.linspace(-radius, radius, self.filtOrder+1) / float(radius)) # lanczos
#win = spsig.hamming(self.filtOrder+1) # sinc-hamming
# this should be automated somehow XXX - idfah
if self.filtOrder <= 6:
cutoff = 2*0.570
elif self.filtOrder <= 8:
cutoff = 2*0.676
elif self.filtOrder <= 12:
cutoff = 2*0.781
elif self.filtOrder <= 16:
cutoff = 2*0.836
elif self.filtOrder <= 32:
cutoff = 2*0.918
elif self.filtOrder <= 64:
cutoff = 2*0.959
# need to fix for multiple pool sizes XXX - idfah
cutoff /= float(self.poolSize)
taps = cutoff * np.linspace(-radius, radius, self.filtOrder+1, dtype=self.dtype)
impulseResponse = cutoff * np.sinc(taps) * win
self.filters = []
nReadoutLayers = 1 if self.nHidden is None else 2
for ni, no in self.layerDims[:-nReadoutLayers]:
noEmb = no*(self.filtOrder+1) # no outs after filter embedding
filtMat = np.zeros(noEmb*2, dtype=self.dtype)
filtMat[noEmb-1::-no] = impulseResponse
# filters strided for embedding
sz = filtMat.itemsize
filtMat = npst.as_strided(filtMat, (no,noEmb), strides=(sz,sz))[::-1].T
self.filters.append(filtMat.copy()) Example 40
def rrc(t):
'''
Input: T, evaluation point (seconds)
Output: value of root-raised-cosine at time T
'''
# Delay between two bits
bit_period = 1/BIT_FREQUENCY
# Total amount of bits to transmit
nb_bits = len(LIST_OF_BITS)
# To be returned (sum of contributions)
s = 0.0
# Max value of rrc
m = 4*BETA/np.pi/np.sqrt(bit_period) + (1-BETA)/np.sqrt(bit_period) + sum(abs(2*rootRaisedCosine(i*bit_period)) for i in range(1, TRUNCATION))
if(t < - TRUNCATION * bit_period or t >= (nb_bits + TRUNCATION) * bit_period):
# T out of support
r = 0.0
else:
# Bits that will affect function at time T
relevant_bits = np.zeros(2*TRUNCATION+1)
for i in range(2*TRUNCATION+1):
j = t/bit_period + i - TRUNCATION
j = int(j) if int(j) <= j else int(j) - 1
if(j >= 0 and j < nb_bits):
relevant_bits[i] = -1 if LIST_OF_BITS[j] == '0' else 1
for i in range(2*TRUNCATION+1):
tt = t/bit_period
tt = t - int(tt)*bit_period if int(tt) <= tt else t - (int(tt)-1)*bit_period
if(t == bit_period * (1 / 4 / BETA + (i - TRUNCATION))):
# L'Hospital's rule because of potential discontinuity
s += relevant_bits[i] * BETA / np.pi / np.sqrt(2*bit_period) * 1 / m * \
((np.pi + 2) * np.sin(np.pi/4/BETA) + \
(np.pi - 2) * np.cos(np.pi/4/BETA))
else:
# General case formula
s += relevant_bits[i] * 4*BETA/np.pi/np.sqrt(bit_period) * 1 / m * \
(np.cos((1 + BETA) * np.pi * ((tt / bit_period - (i-TRUNCATION)))) + \
(1 - BETA) * np.pi / 4 / BETA * \
np.sinc((1 - BETA) * (tt / bit_period - (i-TRUNCATION))))/ \
(1 - (4*BETA*(tt / bit_period - (i-TRUNCATION)))**2)
return s
### ### ### ### ### ### ### Example 41
def noise_power(self, freq):
"""Returns a function to calculate the noise PS at the given freq.
"""
z = freq_to_z(freq)
beam_size = self.beam_size(freq)
A_pix = beam_size**2
A_survey = self.f_sky * 4 * np.pi
tau = A_pix / A_survey * units.year * self.num_year * self.beam_num
# Calculate the comoving size of a frequency bin (at the given freq)
d = self.proper_distance
dxf = (d(freq_to_z(freq - self.freq_width)) - d(z))
# Define the window function in k-space for the parallel and perpendicular directions.
# Use a sinc function for parallel as it is the FT of a top-hat bin. This is probably a bad choice.
def window_par(kpar):
y = kpar * dxf / (4 * np.pi)
return np.sinc(y) * (np.abs(y) < 1.0)
# Azimuthally average over the X and Y window functions to produce an
# overall k_perp window function. Do this by averaging for a set number
# of points k values and then generating an interpolating function to
# appeoximate the full result.
def _int(phi, k):
# Integrand to average over
x = (3e2 * k * d(z)) / (freq * 2 * np.pi)
xx = x * np.cos(phi)
xy = x * np.sin(phi)
return (self.window_x(xx) * self.window_y(xy))**2
def _w_xy_average(k):
# Full averaged window function
return scipy.integrate.fixed_quad(_int, 0, 2 * np.pi, args=(k,), n=1024)[0]
# Generate a log interpolated approximation
k_val = np.linspace(0, self.kmax, 256)
int_val = np.array([_w_xy_average(k)**0.5 for k in k_val])
_w_perp_interp = scipy.interpolate.interp1d(k_val, np.log(int_val))
def window_perp(kperp):
return np.exp(_w_perp_interp(kperp))
# Calculate the comoving volume of a single pixel (beam)
V_pix = A_pix * d(z)**2 * dxf
# Receiver temperature contribution to instrumental Stokes I
T_recv_I = self.T_recv / 2**0.5
return inv_noise_ps_21cm(T_recv_I + self.T_sky(freq), tau, V_pix,
self.freq_width, window_par, window_perp) Example 42
def focus_multiprocessing(self, row):
"""
Focus SAR image with TDBP algorithm. NOTE: Image must be range compressed.
It uses local squint angle (antenna coordinate system) and distances to target to focus.
Parameters
----------
row: int.
image row to be focus.
Returns
-------
list of numpy complex.
List containing entire focused row (numpy complex data) calculated in parallel mode.
"""
# Light speed.
c = 300000000.0
# SAR bandwidth, central frequency and lambda.
sar_B = self.param.get_float_parameter("Radar/B")
sar_f0 = self.param.get_float_parameter("Radar/f0")
sar_lambda = c/sar_f0
nt_fast_time = self.simulated_image.traj.nt
nt_slow_time = self.simulated_image.Nt
# Partial row calculated in parallel mode focusing.
partial_row = np.empty(self.ny, dtype=np.complex128)
x_foc_ind = row
for y_foc_ind in range(self.ny):
foc_lin_ind = x_foc_ind*self.ny + y_foc_ind
# Synthetic range compressed data (matched 2D filter).
# Antenna Enclosure (lobe).
doppler_amplitude_lin = (np.sinc(self.local_squint_ref_traj[foc_lin_ind, :]/self.radar_beamwidth*0.886 ))**2
doppler_amplitude = np.tile(doppler_amplitude_lin, [self.simulated_image.Nt, 1])
# Range amplitude: range positions in raw data of backscattered signal. These are the sincs with range
# migration (range compressed image).
range_amplitude = np.sinc( sar_B*( (np.tile(self.simulated_image.t_axis_fast_time, [nt_fast_time, 1])).transpose()
- np.tile(2*self.distances_ref_traj[foc_lin_ind, :]/c, [nt_slow_time, 1]) ) )
# Limit bandwidth to threshold given by a window. Use only 3dB of antenna lobe for azimuth, limited by squint threshold.
doppler_threshold_win = np.absolute( np.tile(self.local_squint_ref_traj[foc_lin_ind, :], [nt_slow_time, 1]) ) < self.squint_threshold
raw_amplitude = doppler_amplitude*range_amplitude*doppler_threshold_win
# Phase of backscattered signal (2*pi*2*r/lambda).
raw_phase = np.exp(-1j*4*np.pi/sar_lambda*np.tile(self.distances_ref_traj[foc_lin_ind, :], [nt_slow_time, 1]))
# Get module of raw_amplitude (for every xn, yn).
mod_raw_amplitude = np.sum(abs(raw_amplitude)**2)
# Repeat over x,y (slow time and fast time) to normalize.
mod_raw_amplitude = np.tile(mod_raw_amplitude, [nt_slow_time, nt_fast_time])
# Get raw odographer with raw_amplitude and raw_phase, i.e. with amplitude and phase information, and normalize.
raw_to_foc = (np.conjugate(raw_phase))*raw_amplitude/mod_raw_amplitude
partial_row[y_foc_ind] = np.sum(self.new_raw*raw_to_foc)
return list(partial_row) Example 43
def generate_img(self, param):
"""
Generate range compressed image.
Parameters
----------
param: object (ConfigurationManager instance).
ConfigurationManager instance to read parameters from file.
Returns
-------
-.
"""
# Light speed.
c = 300000000.0
# Load beamwidth, bandwidth and central frequency to use locally.
sar_bmw = (param.get_float_parameter("Radar/beamwidth")*np.pi)/180.
sar_B = param.get_float_parameter("Radar/B")
sar_f0 = param.get_float_parameter("Radar/f0")
# Get angles squint and look of view with respect to the antenna coordinate system.
#self.get_angles_antenna()
self.local_look, self.local_squint = Utils.get_angles_antenna(self.traj, self.nom_target)
# Set fast time axis.
start = 2*(min(self.distances))/c - self.fast_time_pixel_margin_mono*self.radar_dt
end = 2*(max(self.distances))/c + self.fast_time_pixel_margin_mono*self.radar_dt
step = self.radar_dt
self.t_axis_fast_time = np.arange(start, end, step)
# Number of elements in fast time axis.
self.Nt = np.size(self.t_axis_fast_time)
self.freq_axis_fftshift = Utils.freq_axis(self.radar_dt, self.Nt, False, True)
sar_lambda = c/sar_f0
# Doppler amplitude (envolvente de la antena).
doppler_amplitude = (np.sinc( (np.tile(self.local_squint, [self.Nt, 1]))/sar_bmw*(2*0.443) ))**2
# Range amplitude: range positions in raw data of backscattered signal.
Nd = np.size(self.distances)
range_amplitude = np.sinc( sar_B*( (np.tile(self.t_axis_fast_time, [Nd, 1])).transpose() - np.tile(2*self.distances/c, [self.Nt, 1]) ) )
# Signal phase received: 2*pi*2*r/lambda.
signal_phase = np.exp(-1j*4*np.pi/sar_lambda*np.tile(self.distances, [self.Nt, 1]))
# Generate range compressed simulated image.
self.image = doppler_amplitude*range_amplitude*signal_phase 