def _lagged_coherence_1freq(x, f, Fs, N_cycles=3, f_step=1):
"""Calculate lagged coherence of x at frequency f using the hanning-taper FFT method"""
Nsamp = int(np.ceil(N_cycles*Fs / f))
# For each N-cycle chunk, calculate phase
chunks = _nonoverlapping_chunks(x,Nsamp)
C = len(chunks)
hann_window = signal.hanning(Nsamp)
fourier_f = np.fft.fftfreq(Nsamp,1/float(Fs))
fourier_f_idx = _arg_closest_value(fourier_f,f)
fourier_coefsoi = np.zeros(C,dtype=complex)
for i2, c in enumerate(chunks):
fourier_coef = np.fft.fft(c*hann_window)
fourier_coefsoi[i2] = fourier_coef[fourier_f_idx]
lcs_num = 0
for i2 in range(C-1):
lcs_num += fourier_coefsoi[i2]*np.conj(fourier_coefsoi[i2+1])
lcs_denom = np.sqrt(np.sum(np.abs(fourier_coefsoi[:-1])**2)*np.sum(np.abs(fourier_coefsoi[1:])**2))
return np.abs(lcs_num/lcs_denom)
python类conj()的实例源码
def _create_rotational_weights_for_elements(self, kpoint, transformation_matrix, vectors):
"""
Parameters
----------
kpoint : 1d array
Reciprocal space point in fractional coordinates for PC.
vectors : (..., natoms_p * ndims, nbands) array
Vectors for SC after translational projection.
"""
projected_vectors = self._rotational_projector.project_vectors(
vectors, kpoint, transformation_matrix)
nirreps, natoms_p, nelms, tmp, nbands = projected_vectors.shape
shape = (nirreps, natoms_p, nelms, natoms_p, nelms, nbands)
weights = np.zeros(shape, dtype=complex)
for i in range(nirreps):
for j in range(nbands):
weights[i, ..., j] = np.inner(
np.conj(projected_vectors[i, ..., j]), projected_vectors[i, ..., j])
return weights, projected_vectors
def get_power_spectral_density_matrix(observation, mask=None):
"""
Calculates the weighted power spectral density matrix.
This does not yet work with more than one target mask.
:param observation: Complex observations with shape (bins, sensors, frames)
:param mask: Masks with shape (bins, frames) or (bins, 1, frames)
:return: PSD matrix with shape (bins, sensors, sensors)
"""
bins, sensors, frames = observation.shape
if mask is None:
mask = np.ones((bins, frames))
if mask.ndim == 2:
mask = mask[:, np.newaxis, :]
normalization = np.maximum(np.sum(mask, axis=-1, keepdims=True), 1e-6)
psd = np.einsum('...dt,...et->...de', mask * observation,
observation.conj())
psd /= normalization
return psd
def get_mvdr_vector(atf_vector, noise_psd_matrix):
"""
Returns the MVDR beamforming vector.
:param atf_vector: Acoustic transfer function vector
with shape (..., bins, sensors)
:param noise_psd_matrix: Noise PSD matrix
with shape (bins, sensors, sensors)
:return: Set of beamforming vectors with shape (..., bins, sensors)
"""
while atf_vector.ndim > noise_psd_matrix.ndim - 1:
noise_psd_matrix = np.expand_dims(noise_psd_matrix, axis=0)
# Make sure matrix is hermitian
noise_psd_matrix = 0.5 * (
noise_psd_matrix + np.conj(noise_psd_matrix.swapaxes(-1, -2)))
numerator = solve(noise_psd_matrix, atf_vector)
denominator = np.einsum('...d,...d->...', atf_vector.conj(), numerator)
beamforming_vector = numerator / np.expand_dims(denominator, axis=-1)
return beamforming_vector
def autocorrelation(self):
"Autocorrelation as a function of time"
if self.__autocorrelation is not None:
return self.__autocorrelationTimeSeries, self.__autocorrelation
negT = -np.flipud(self.timeSeries[1:])
autocorrelationTime = np.hstack((negT, self.timeSeries))
self.__autocorrelationTimeSeries = autocorrelationTime
initialWF = self[0]
ACF = []
for WF in self:
ACF.append(WF.overlap(initialWF))
ACF = np.array(ACF)
negACF = np.conj(np.flipud(ACF[1:]))
totalACF = np.hstack((negACF, ACF))
self.__autocorrelation = totalACF
return self.__autocorrelationTimeSeries, self.__autocorrelation
def swap_Nq(fft_y, fu, fft_x, N):
f = fu[:, 0].copy()
fft_x[0] = f[0].real
fft_x[1:N//2] = 0.5*(f[1:N//2] + np.conj(f[:N//2:-1]))
fft_x[N//2] = f[N//2].real
fu[:N//2+1, 0] = fft_x[:N//2+1]
fu[N//2+1:, 0] = np.conj(fft_x[(N//2-1):0:-1])
fft_y[0] = f[0].imag
fft_y[1:N//2] = -0.5*1j*(f[1:N//2] - np.conj(f[:N//2:-1]))
fft_y[N//2] = f[N//2].imag
fft_y[N//2+1:] = np.conj(fft_y[(N//2-1):0:-1])
return fft_y
def nufft_alpha_kb_fit(N, J, K):
'''
find out parameters alpha and beta
of scaling factor st['sn']
Note, when J = 1 , alpha is hardwired as [1,0,0...]
(uniform scaling factor)
'''
beta = 1
Nmid = (N - 1.0) / 2.0
if N > 40:
L = 13
else:
L = numpy.ceil(N / 3)
nlist = numpy.arange(0, N) * 1.0 - Nmid
(kb_a, kb_m) = kaiser_bessel('string', J, 'best', 0, K / N)
if J > 1:
sn_kaiser = 1 / kaiser_bessel_ft(nlist / K, J, kb_a, kb_m, 1.0)
elif J == 1: # The case when samples are on regular grids
sn_kaiser = numpy.ones((1, N), dtype=dtype)
gam = 2 * numpy.pi / K
X_ant = beta * gam * nlist.reshape((N, 1), order='F')
X_post = numpy.arange(0, L + 1)
X_post = X_post.reshape((1, L + 1), order='F')
X = numpy.dot(X_ant, X_post) # [N,L]
X = numpy.cos(X)
sn_kaiser = sn_kaiser.reshape((N, 1), order='F').conj()
X = numpy.array(X, dtype=dtype)
sn_kaiser = numpy.array(sn_kaiser, dtype=dtype)
coef = numpy.linalg.lstsq(numpy.nan_to_num(X), numpy.nan_to_num(sn_kaiser))[0] # (X \ sn_kaiser.H);
alphas = coef
if J > 1:
alphas[0] = alphas[0]
alphas[1:] = alphas[1:] / 2.0
elif J == 1: # cases on grids
alphas[0] = 1.0
alphas[1:] = 0.0
alphas = numpy.real(alphas)
return (alphas, beta)
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 series_nfft(series,
oversample=4):
"""
note that output period units are [days] (so is frequency)
"""
M = len(series)
if not is_power_of_2(M):
raise ValueError('series length {} is not a power of 2'.format(len(series)))
N = M
if N % 2 == 1:
# number of frequency samples must be even
N += 1
N *= oversample
# re-grid time the interval [-1/2, 1/2) as required by nfft
time = series.index.astype(NP.int64) / 1e9
time -= time[0]
b = -0.5
a = (M - 1) / (M * time[-1])
x = a * time + b
# setup for nfft computation
plan = NFFT(N, M)
plan.x = x
plan.f = series.values
plan.precompute()
# compute nfft (note that conjugation is necessary because of the
# difference in transform sign convention)
x_nfft = NP.conj(plan.adjoint())
# calculate frequencies and periods
dt = ((series.index[-1] - series.index[0]) / M).total_seconds() / DAYS_TO_SECONDS
f_range = NP.fft.fftshift(NP.fft.fftfreq(N, dt))
T_range = 1 / f_range
return x_nfft, f_range, T_range
def spatFT(data, position_grid, order_max=10, spherical_harmonic_bases=None):
''' Spatial Fourier Transform
Parameters
----------
data : array_like
Data to be transformed, with signals in rows and frequency bins in columns
order_max : int, optional
Maximum transform order (Default: 10)
position_grid : array_like or io.SphericalGrid
Azimuths/Colatitudes/Gridweights of spatial sampling points
Returns
-------
Pnm : array_like
Spatial Fourier Coefficients with nm coeffs in rows and FFT bins in columns
'''
data = _np.atleast_2d(data)
number_of_signals, FFTLength = data.shape
position_grid = SphericalGrid(*position_grid)
# Re-generate spherical harmonic bases if they were not provided or their order is too small
if (spherical_harmonic_bases is None or
spherical_harmonic_bases.shape[0] < number_of_signals or
spherical_harmonic_bases.shape[1] < (order_max + 1) ** 2):
spherical_harmonic_bases = sph_harm_all(order_max, position_grid.azimuth, position_grid.colatitude)
spherical_harmonic_bases = (_np.conj(spherical_harmonic_bases).T * (4 * _np.pi * position_grid.weight))
return _np.inner(spherical_harmonic_bases, data.T)
def _random_state(sites, ldim, randstate=None):
"""Returns a random positive semidefinite operator of shape (ldim, ldim) *
sites normalized to Tr rho = 1, i.e. a mixed state with local dimension
`ldim` living on `sites` sites. Note that the returned state is positive
semidefinite only when interpreted in global form (see
:func:`tools.global_to_local`)
:param sites: Number of local sites
:param ldim: Local ldimension
:param randstate: numpy.random.RandomState instance or None
:returns: numpy.ndarray of shape (ldim, ldim) * sites
>>> from numpy.linalg import eigvalsh
>>> rho = _random_state(3, 2).reshape((2**3, 2**3))
>>> all(eigvalsh(rho) >= 0)
True
>>> np.abs(np.trace(rho) - 1) < 1e-6
True
"""
shape = (ldim**sites, ldim**sites)
mat = _zrandn(shape, randstate=randstate)
rho = np.conj(mat.T).dot(mat)
rho /= np.trace(rho)
return rho.reshape((ldim,) * 2 * sites)
####################################
# Factory functions for MPArrays #
####################################
def random_mpo(sites, ldim, rank, randstate=None, hermitian=False,
normalized=True, force_rank=False):
"""Returns an hermitian MPO with randomly choosen local tensors
:param sites: Number of sites
:param ldim: Local dimension
:param rank: Rank
:param randstate: numpy.random.RandomState instance or None
:param hermitian: Is the operator supposed to be hermitian
:param normalized: Operator should have unit norm
:param force_rank: If True, the rank is exaclty `rank`.
Otherwise, it might be reduced if we reach the maximum sensible rank.
:returns: randomly choosen matrix product operator
>>> mpo = random_mpo(4, 2, 10, force_rank=True)
>>> mpo.ranks, mpo.shape
((10, 10, 10), ((2, 2), (2, 2), (2, 2), (2, 2)))
>>> mpo.canonical_form
(0, 4)
"""
mpo = random_mpa(sites, (ldim,) * 2, rank, randstate=randstate,
force_rank=force_rank, dtype=np.complex_)
if hermitian:
# make mpa Herimitan in place, without increasing rank:
ltens = (l + l.swapaxes(1, 2).conj() for l in mpo.lt)
mpo = mp.MPArray(ltens)
if normalized:
# we do this with a copy to ensure the returned state is not
# normalized
mpo /= mp.norm(mpo.copy())
return mpo
def random_mpdo(sites, ldim, rank, randstate=np.random):
"""Returns a randomly choosen matrix product density operator (i.e.
positive semidefinite matrix product operator with trace 1).
:param sites: Number of sites
:param ldim: Local dimension
:param rank: Rank
:param randstate: numpy.random.RandomState instance
:returns: randomly choosen classicaly correlated matrix product density op.
>>> rho = random_mpdo(4, 2, 4)
>>> rho.ranks, rho.shape
((4, 4, 4), ((2, 2), (2, 2), (2, 2), (2, 2)))
>>> rho.canonical_form
(0, 4)
"""
# generate density matrix as a mixture of `rank` pure product states
psis = [random_mps(sites, ldim, 1, randstate=randstate) for _ in range(rank)]
weights = (lambda x: x / np.sum(x))(randstate.rand(rank))
rho = mp.sumup(mpsmpo.mps_to_mpo(psi) * weight
for weight, psi in zip(weights, psis))
# Scramble the local tensors
for n, rank in enumerate(rho.ranks):
unitary = _unitary_haar(rank, randstate)
rho.lt[n] = matdot(rho.lt[n], unitary)
rho.lt[n + 1] = matdot(np.transpose(unitary).conj(), rho.lt[n + 1])
rho /= mp.trace(rho)
return rho
def test_conjugations(nr_sites, local_dim, _, rgen, dtype):
op = factory._random_op(nr_sites, local_dim, randstate=rgen, dtype=dtype)
mpo = mp.MPArray.from_array(op, 2)
assert_array_almost_equal(np.conj(op), mpo.conj().to_array())
assert mpo.conj().dtype == dtype
mpo.canonicalize()
mpo_c = mpo.conj()
assert_correct_normalization(mpo_c)
def test_norm(nr_sites, local_dim, rank, dtype, rgen):
mp_psi = factory.random_mpa(nr_sites, local_dim, rank, randstate=rgen,
dtype=dtype)
psi = mp_psi.to_array()
assert_almost_equal(mp.inner(mp_psi, mp_psi), mp.norm(mp_psi)**2)
assert_almost_equal(np.sum(psi.conj() * psi), mp.norm(mp_psi)**2)
def SLshearrecadjoint2D(X, shearletSystem):
"""
Adjoint of (pseudo-)inverse of 2D data.
Note that this is also the (pseudo-)inverse of the adjoint.
Usage:
coeffs = SLshearrecadjoint2D(X, shearletSystem)
Input:
X : 2D data.
shearletSystem: Structure containing a shearlet system. This
should be the same system as the one
previously used for decomposition.
Output:
coeffs: X x Y x N array of shearlet coefficients.
"""
# skipping useGPU stuff...
# STUFF
Xfreq = np.divide(fftlib.fftshift(fftlib.fft2(fftlib.ifftshift(X))), shearletSystem["dualFrameWeights"])
coeffs = np.zeros(shearletSystem["shearlets"].shape, dtype=complex)
for j in range(shearletSystem["nShearlets"]):
coeffs[:,:,j] = fftlib.fftshift(fftlib.ifft2(fftlib.ifftshift(Xfreq*np.conj(shearletSystem["shearlets"][:,:,j]))))
return np.real(coeffs).astype(X.dtype)
#
##############################################################################
def test_simple_conjugate(self):
ref = np.conj(np.sqrt(np.complex(1, 1)))
def f(z):
return np.sqrt(np.conj(z))
yield check_complex_value, f, 1, 1, ref.real, ref.imag, False
#def test_branch_cut(self):
# _check_branch_cut(f, -1, 0, 1, -1)
def test_cabs_inf_nan(self):
x, y = [], []
# cabs(+-nan + nani) returns nan
x.append(np.nan)
y.append(np.nan)
yield check_real_value, np.abs, np.nan, np.nan, np.nan
x.append(np.nan)
y.append(-np.nan)
yield check_real_value, np.abs, -np.nan, np.nan, np.nan
# According to C99 standard, if exactly one of the real/part is inf and
# the other nan, then cabs should return inf
x.append(np.inf)
y.append(np.nan)
yield check_real_value, np.abs, np.inf, np.nan, np.inf
x.append(-np.inf)
y.append(np.nan)
yield check_real_value, np.abs, -np.inf, np.nan, np.inf
# cabs(conj(z)) == conj(cabs(z)) (= cabs(z))
def f(a):
return np.abs(np.conj(a))
def g(a, b):
return np.abs(np.complex(a, b))
xa = np.array(x, dtype=np.complex)
for i in range(len(xa)):
ref = g(x[i], y[i])
yield check_real_value, f, x[i], y[i], ref
def cor_fft(x1,x2,sigma):
dist11 = np.sum(np.square(x1))
dist22 = np.sum(np.square(x2))
if len(x1.shape)==2:
c = np.fft.ifft2((np.conj(fft(x1))*fft(x2)))
else:
c = np.fft.ifft2(np.sum(np.conj(fft(x1))*fft(x2),2))
dist= dist11-2*c+dist22
cor = np.exp(-1*dist/(sigma**2*x1.size))
cor = np.real(cor)
return cor
