def _ncc_c(x, y):
"""
>>> _ncc_c([1,2,3,4], [1,2,3,4])
array([ 0.13333333, 0.36666667, 0.66666667, 1. , 0.66666667,
0.36666667, 0.13333333])
>>> _ncc_c([1,1,1], [1,1,1])
array([ 0.33333333, 0.66666667, 1. , 0.66666667, 0.33333333])
>>> _ncc_c([1,2,3], [-1,-1,-1])
array([-0.15430335, -0.46291005, -0.9258201 , -0.77151675, -0.46291005])
"""
den = np.array(norm(x) * norm(y))
den[den == 0] = np.Inf
x_len = len(x)
fft_size = 1<<(2*x_len-1).bit_length()
cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size)))
cc = np.concatenate((cc[-(x_len-1):], cc[:x_len]))
return np.real(cc) / den
python类ifft()的实例源码
def cconv(a, b):
"""
Circular convolution of vectors
Computes the circular convolution of two vectors a and b via their
fast fourier transforms
a \ast b = \mathcal{F}^{-1}(\mathcal{F}(a) \odot \mathcal{F}(b))
Parameter
---------
a: real valued array (shape N)
b: real valued array (shape N)
Returns
-------
c: real valued array (shape N), representing the circular
convolution of a and b
"""
return ifft(fft(a) * fft(b)).real
def ccorr(a, b):
"""
Circular correlation of vectors
Computes the circular correlation of two vectors a and b via their
fast fourier transforms
a \ast b = \mathcal{F}^{-1}(\overline{\mathcal{F}(a)} \odot \mathcal{F}(b))
Parameter
---------
a: real valued array (shape N)
b: real valued array (shape N)
Returns
-------
c: real valued array (shape N), representing the circular
correlation of a and b
"""
return ifft(np.conj(fft(a)) * fft(b)).real
def hilbert_fredriksen(f, ker=None):
"""
Performs Hilbert transform of f.
Parameters
----------
f : array
Values of function on a equidistant grid.
ker : array
Kernel used when performing Hilbert transform using FFT.
Returns
-------
array
Hilbert transform of f.
"""
if ker is None:
ker = kernel_fredriksen(len(f))
n = len(f)
fpad = fft(np.concatenate( (f,np.zeros(len(ker)-n)) ))
r = ifft(fpad*ker)
return r[0:n]
def cross_correlation_using_fft(self, x, y):
f1 = fft(x)
f2 = fft(np.flipud(y))
cc = np.real(ifft(f1 * f2))
return fftshift(cc)
# clean up
# def close(self):
# close serial
# self.ser.flush()
# self.ser.close()
# Main
def init_pyfftw(x, effort=effort[0], wis=False):
N = len(x)
a = pyfftw.n_byte_align_empty(int(N), n, 'complex128')
a[:] = x
if wis is not False:
pyfftw.import_wisdom(wis)
fft = pyfftw.builders.fft(a, threads=8)
ifft = pyfftw.builders.ifft(a, threads=8)
else:
fft = pyfftw.builders.fft(a, planner_effort=effort, threads=8)
ifft = pyfftw.builders.ifft(a, planner_effort=effort, threads=8)
return fft, ifft
def delayseq(x, delay_sec:float, fs:int):
"""
x: input 1-D signal
delay_sec: amount to shift signal [seconds]
fs: sampling frequency [Hz]
xs: time-shifted signal
"""
assert x.ndim == 1, 'only 1-D signals for now'
delay_samples = delay_sec*fs
delay_int = round(delay_samples)
nfft = nextpow2(x.size+delay_int)
fbins = 2*pi*ifftshift((arange(nfft)-nfft//2))/nfft
X = fft(x,nfft)
Xs = ifft(X*exp(-1j*delay_samples*fbins))
if isreal(x[0]):
Xs = Xs.real
xs = zeros_like(x)
xs[delay_int:] = Xs[delay_int:x.size]
return xs
def invFourier(X,fs,N): # Approximate inverse Fourier transform on the interval (-T/2,T/2)
NFFT=len(X)
k=np.arange(NFFT)
x=fs*np.exp(1j*np.pi*(N/2.-k))*ifft(X*np.exp(-1j*np.pi*k*N/NFFT),NFFT)
x=x[:N]
return x
