我们从Python开源项目中,提取了以下43个代码示例,用于说明如何使用numpy.sinc()。
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
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 #----------------------------------
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
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 : ####################################################################
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
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
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
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
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]
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
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
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
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
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)
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)
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)))
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):
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
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
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
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
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
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
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
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
def lanczos(n, radius=3): taps = np.linspace(-radius,radius, n) return np.sinc(taps/radius)
def sinc(n, radius=3, freq=1.0): taps = np.linspace(-radius,radius, n) return freq * np.sinc(freq*taps)
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)
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) + '.')
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))
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
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)
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
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)
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
def dirichlet(x): return numpy.sinc(x)
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())
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 ### ### ### ### ### ### ###
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)
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)
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