# Python numpy.fft.irfft() Examples

The following are 25 code examples for showing how to use numpy.fft.irfft(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module , or try the search function .

Example 1
def pink(N):
"""
Pink noise.

:param N: Amount of samples.

Pink noise has equal power in bands that are proportionally wide.
Power density decreases with 3 dB per octave.

"""
# This method uses the filter with the following coefficients.
# b = np.array([0.049922035, -0.095993537, 0.050612699, -0.004408786])
# a = np.array([1, -2.494956002, 2.017265875, -0.522189400])
# return lfilter(B, A, np.random.randn(N))
# Another way would be using the FFT
# x = np.random.randn(N)
# X = rfft(x) / N
uneven = N % 2
X = np.random.randn(N // 2 + 1 + uneven) + 1j * np.random.randn(N // 2 + 1 + uneven)
S = np.sqrt(np.arange(len(X)) + 1.)  # +1 to avoid divide by zero
y = (irfft(X / S)).real
if uneven:
y = y[:-1]
return normalize(y) 
Example 2
def mass(query, ts):
"""
>>> np.round(mass(np.array([0.0, 1.0, -1.0, 0.0]), np.array([-1, 1, 0, 0, -1, 1])), 3)
array([ 2.   ,  2.828,  2.   ])
"""
query = zNormalize(query)
m = len(query)
n = len(ts)

stdv = movstd(ts, m)
query = query[::-1]
query = np.pad(query, (0, n - m), 'constant')
dots = fft.irfft(fft.rfft(ts) * fft.rfft(query))
return np.sqrt(2 * (m - (dots[m - 1 :] / stdv))) 
Example 3
def impedance2wake(f, y):
"""
Fourier transform with exp(-iwt)
f -    Hz
y -    Om
s    - Meter
w -    V/C
"""
df = f[1] - f[0]
n = len(f)
s = 1./df*np.arange(n)/n*speed_of_light
w = n*df*irfft(y, n)
return s, w 
Example 4
def bench_random(self):
from numpy.fft import irfft as numpy_irfft

print()
print('Inverse Fast Fourier Transform (real data)')
print('==================================')
print(' size |  scipy  |  numpy  ')
print('----------------------------------')
for size,repeat in [(100,7000),(1000,2000),
(256,10000),
(512,10000),
(1024,1000),
(2048,1000),
(2048*2,500),
(2048*4,500),
]:
print('%5s' % size, end=' ')
sys.stdout.flush()

x = random([size]).astype(double)
x1 = zeros(size/2+1,dtype=cdouble)
x1[0] = x[0]
for i in range(1,size/2):
x1[i] = x[2*i-1] + 1j * x[2*i]
if not size % 2:
x1[-1] = x[-1]
y = irfft(x)

print('|%8.2f' % measure('irfft(x)',repeat), end=' ')
sys.stdout.flush()

assert_array_almost_equal(numpy_irfft(x1,size),y)
print('|%8.2f' % measure('numpy_irfft(x1,size)',repeat), end=' ')
sys.stdout.flush()

print(' (secs for %s calls)' % (repeat))

sys.stdout.flush() 
Example 5
def test_fourier_gaussian_real01(self):
for shape in [(32, 16), (31, 15)]:
for type in [numpy.float32, numpy.float64]:
a = numpy.zeros(shape, type)
a[0, 0] = 1.0
a = fft.rfft(a, shape[0], 0)
a = fft.fft(a, shape[1], 1)
a = ndimage.fourier_gaussian(a, [5.0, 2.5],
shape[0], 0)
a = fft.ifft(a, shape[1], 1)
a = fft.irfft(a, shape[0], 0)
assert_almost_equal(ndimage.sum(a), 1) 
Example 6
def test_fourier_uniform_real01(self):
for shape in [(32, 16), (31, 15)]:
for type in [numpy.float32, numpy.float64]:
a = numpy.zeros(shape, type)
a[0, 0] = 1.0
a = fft.rfft(a, shape[0], 0)
a = fft.fft(a, shape[1], 1)
a = ndimage.fourier_uniform(a, [5.0, 2.5],
shape[0], 0)
a = fft.ifft(a, shape[1], 1)
a = fft.irfft(a, shape[0], 0)
assert_almost_equal(ndimage.sum(a), 1.0) 
Example 7
def test_fourier_shift_real01(self):
for shape in [(32, 16), (31, 15)]:
for dtype in [numpy.float32, numpy.float64]:
expected = numpy.arange(shape[0] * shape[1], dtype=dtype)
expected.shape = shape
a = fft.rfft(expected, shape[0], 0)
a = fft.fft(a, shape[1], 1)
a = ndimage.fourier_shift(a, [1, 1], shape[0], 0)
a = fft.ifft(a, shape[1], 1)
a = fft.irfft(a, shape[0], 0)
assert_array_almost_equal(a[1:, 1:], expected[:-1, :-1])
assert_array_almost_equal(a.imag, numpy.zeros(shape)) 
Example 8
def test_fourier_ellipsoid_real01(self):
for shape in [(32, 16), (31, 15)]:
for type in [numpy.float32, numpy.float64]:
a = numpy.zeros(shape, type)
a[0, 0] = 1.0
a = fft.rfft(a, shape[0], 0)
a = fft.fft(a, shape[1], 1)
a = ndimage.fourier_ellipsoid(a, [5.0, 2.5],
shape[0], 0)
a = fft.ifft(a, shape[1], 1)
a = fft.irfft(a, shape[0], 0)
assert_almost_equal(ndimage.sum(a), 1.0) 
Example 9
def test_fourier_gaussian_real01(self):
for shape in [(32, 16), (31, 15)]:
for type_, dec in zip([numpy.float32, numpy.float64], [6, 14]):
a = numpy.zeros(shape, type_)
a[0, 0] = 1.0
a = fft.rfft(a, shape[0], 0)
a = fft.fft(a, shape[1], 1)
a = ndimage.fourier_gaussian(a, [5.0, 2.5], shape[0], 0)
a = fft.ifft(a, shape[1], 1)
a = fft.irfft(a, shape[0], 0)
assert_almost_equal(ndimage.sum(a), 1, decimal=dec) 
Example 10
def test_fourier_uniform_real01(self):
for shape in [(32, 16), (31, 15)]:
for type_, dec in zip([numpy.float32, numpy.float64], [6, 14]):
a = numpy.zeros(shape, type_)
a[0, 0] = 1.0
a = fft.rfft(a, shape[0], 0)
a = fft.fft(a, shape[1], 1)
a = ndimage.fourier_uniform(a, [5.0, 2.5], shape[0], 0)
a = fft.ifft(a, shape[1], 1)
a = fft.irfft(a, shape[0], 0)
assert_almost_equal(ndimage.sum(a), 1.0, decimal=dec) 
Example 11
def test_fourier_shift_real01(self):
for shape in [(32, 16), (31, 15)]:
for type_, dec in zip([numpy.float32, numpy.float64], [4, 11]):
expected = numpy.arange(shape[0] * shape[1], dtype=type_)
expected.shape = shape
a = fft.rfft(expected, shape[0], 0)
a = fft.fft(a, shape[1], 1)
a = ndimage.fourier_shift(a, [1, 1], shape[0], 0)
a = fft.ifft(a, shape[1], 1)
a = fft.irfft(a, shape[0], 0)
assert_array_almost_equal(a[1:, 1:], expected[:-1, :-1],
decimal=dec)
assert_array_almost_equal(a.imag, numpy.zeros(shape),
decimal=dec) 
Example 12
def test_fourier_ellipsoid_real01(self):
for shape in [(32, 16), (31, 15)]:
for type_, dec in zip([numpy.float32, numpy.float64], [5, 14]):
a = numpy.zeros(shape, type_)
a[0, 0] = 1.0
a = fft.rfft(a, shape[0], 0)
a = fft.fft(a, shape[1], 1)
a = ndimage.fourier_ellipsoid(a, [5.0, 2.5],
shape[0], 0)
a = fft.ifft(a, shape[1], 1)
a = fft.irfft(a, shape[0], 0)
assert_almost_equal(ndimage.sum(a), 1.0, decimal=dec) 
Example 13
def lpc_python(x, order):
def levinson(R,order):
""" input: autocorrelation and order, output: LPC coefficients """
a   = np.zeros(order+2)
a[0] = -1
k = np.zeros(order+1)
# step 1: initialize prediction error "e" to R[0]
e = R[0]
# step 2: iterate over [1:order]
for i in range(1,order+1):
# step 2-1: calculate PARCOR coefficients
k[i]= (R[i] - np.sum(a[1:i] * R[i-1:0:-1])) / e
# step 2-2: update LPCs
a[i] = np.copy(k[i])
a_old = np.copy(a)
for j in range(1,i):
a[j] -=  k[i]* a_old[i-j]
# step 2-3: update prediction error "e"
e = e * (1.0 - k[i]**2)
return -1*a[0:order+1], e, -1*k[1:]
# N: compute next power of 2
n = x.shape[0]
N = int(np.power(2, np.ceil(np.log2(2 * n - 1))))
# calculate autocorrelation using FFT
X = rfft(x, N)
r = irfft(abs(X) ** 2)[:order+1]
return levinson(r, order) 
Example 14
def multS(s, v, L, TP=False):

N = s.shape[0]
vp = prepare_v(v, N, L, TP=TP)
p = irfft(rfft(vp) * rfft(s))
if not TP:
return p[:L]
return p[L - 1:] 
Example 15
def autocorr(x):
""" Fast autocorrelation computation using the FFT """

X = fft.rfft(x, n=2*x.shape[0])
r = fft.irfft(np.abs(X)**2)

return r[:x.shape[0]] / (x.shape[0] - 1) 
Example 16
def toeplitz_multiplication(c, r, A, **kwargs):
"""
Compute numpy.dot(scipy.linalg.toeplitz(c,r), A) using the FFT.

Parameters
----------
c: ndarray
the first column of the Toeplitz matrix
r: ndarray
the first row of the Toeplitz matrix
A: ndarray
the matrix to multiply on the right
"""

m = c.shape[0]
n = r.shape[0]

fft_len = int(2**np.ceil(np.log2(m+n-1)))
zp = fft_len - m - n + 1

if A.shape[0] != n:
raise ValueError('A dimensions not compatible with toeplitz(c,r)')

x = np.concatenate((c, np.zeros(zp, dtype=c.dtype), r[-1:0:-1]))
xf = np.fft.rfft(x, n=fft_len)

Af = np.fft.rfft(A, n=fft_len, axis=0)

return np.fft.irfft((Af.T*xf).T, n=fft_len, axis=0)[:m,] 
Example 17
def synthesis(self, X=None):
"""
Perform time synthesis of frequency domain to real signal using the
inverse DFT.

Parameters
----------
X : numpy array, optional
Complex signal in frequency domain. Default is to use DFT computed
from analysis.

Returns
-------
numpy array
IDFT of self.X or input if given.
"""

# check for valid input
if X is not None:
if X.shape != self.X.shape:
raise ValueError('Invalid input dimensions! Got (%d, %d), expecting (%d, %d).'
% (X.shape[0], X.shape[1],
self.X.shape[0], self.X.shape[1]))

self.X[:,] = X

# inverse DFT
if self.transform == 'fftw':
self.b[:] = self.X
self.x[:,] = self._backward()
elif self.transform == 'mkl':
self.x[:,] = mkl_fft.irfft_numpy(self.X, self.nfft, axis=self.axis)
else:
self.x[:,] = irfft(self.X, self.nfft, axis=self.axis)

# apply window if needed
if self.synthesis_window is not None:
np.multiply(self.synthesis_window, self.x, self.x)

return self.x 
Example 18
def deconvolve(y, s, length=None, thresh=0.0):
'''
Deconvolve an excitation signal from an impulse response

Parameters
----------

y : ndarray
The recording
s : ndarray
The excitation signal
length: int, optional
the length of the impulse response to deconvolve
thresh : float, optional
ignore frequency bins with power lower than this
'''

# FFT length including zero padding
n = y.shape[0] + s.shape[0] - 1

# let FFT size be even for convenience
if n % 2 != 0:
n += 1

# when unknown, pick the filter size as size of test signal
if length is None:
length = n

# Forward transforms
Y  = fft.rfft(np.array(y, dtype=np.float32), n=n) / np.sqrt(n)
S = fft.rfft(np.array(s, dtype=np.float32), n=n) / np.sqrt(n)

# Only do the division where S is large enough
H = np.zeros(*Y.shape, dtype=Y.dtype)
I = np.where(np.abs(S) > thresh)
H[I] = Y[I] / S[I]

# Inverse transform
h = fft.irfft(H, n=n)

return h[:length] 
Example 19
def correlate(x1, x2, interp=1, phat=False):
'''
Compute the cross-correlation between x1 and x2

Parameters
----------
x1,x2: array_like
The data arrays
interp: int, optional
The interpolation factor for the output array, default 1.
phat: bool, optional
Apply the PHAT weighting (default False)

Returns
-------
The cross-correlation between the two arrays
'''

N1 = x1.shape[0]
N2 = x2.shape[0]

N = N1 + N2 - 1

X1 = fft.rfft(x1, n=N)
X2 = fft.rfft(x2, n=N)

if phat:
eps1 = np.mean(np.abs(X1)) * 1e-10
X1 /= (np.abs(X1) + eps1)
eps2 = np.mean(np.abs(X2)) * 1e-10
X2 /= (np.abs(X2) + eps2)

m = np.minimum(N1, N2)

out = fft.irfft(X1*np.conj(X2), n=int(N*interp))

return np.concatenate([out[-interp*(N2-1):], out[:(interp*N1)]]) 
Example 20
def slidingDotProduct(query,ts):
"""
Calculate the dot product between a query and all subsequences of length(query) in the timeseries ts. Note that we use Numpy's rfft method instead of fft.

Parameters
----------
query: Specific time series query to evaluate.
ts: Time series to calculate the query's sliding dot product against.
"""

m = len(query)
n = len(ts)

#If length is odd, zero-pad time time series
if n%2 ==1:
ts = np.insert(ts,0,0)

#If length is odd, zero-pad query
if m%2 == 1:
query = np.insert(query,0,0)

#This reverses the array
query = query[::-1]

#Determine trim length for dot product. Note that zero-padding of the query has no effect on array length, which is solely determined by the longest vector

dot_product = fft.irfft(fft.rfft(ts)*fft.rfft(query))

#Note that we only care about the dot product results from index m-1 onwards, as the first few values aren't true dot products (due to the way the FFT works for dot products)
return dot_product[trim :] 
Example 21
def generate(self):
X = np.random.randn(self.N // 2 + 1 + self.uneven) + 1j * np.random.randn(self.N // 2 + 1 + self.uneven)
y = (irfft(X / self.S)).real
if self.uneven:
y = y[:-1]
return normalize(y) 
Example 22
def mkl_toeplitz_multiplication(c, r, A, A_padded=False, out=None, fft_len=None):
"""
Compute numpy.dot(scipy.linalg.toeplitz(c,r), A) using the FFT from the mkl_fft package.

Parameters
----------
c: ndarray
the first column of the Toeplitz matrix
r: ndarray
the first row of the Toeplitz matrix
A: ndarray
the matrix to multiply on the right
the A matrix can be pre-padded with zeros by the user, if this is the case
set to True
out: ndarray, optional
an ndarray to store the output of the multiplication
fft_len: int, optional
specify the length of the FFT to use
"""

if not has_mklfft:
raise ValueError('Import mkl_fft package unavailable. Install from https://github.com/LCAV/mkl_fft')

m = c.shape[0]
n = r.shape[0]

if fft_len is None:
fft_len = int(2**np.ceil(np.log2(m+n-1)))
zp = fft_len - m - n + 1

if (not A_padded and A.shape[0] != n) or (A_padded and A.shape[0] != fft_len):
raise ValueError('A dimensions not compatible with toeplitz(c,r)')

x = np.concatenate((c, np.zeros(zp, dtype=c.dtype), r[-1:0:-1]))
xf = fft.rfft(x, n=fft_len)

if out is not None:
fft.rfft(A, n=fft_len, axis=0, out=out)
else:
out = fft.rfft(A, n=fft_len, axis=0)

out *= xf[:,None]

fft.irfft(out, n=fft_len, axis=0, out=A)
else:
A = fft.irfft(out, n=fft_len, axis=0)

return A[:m,] 
Example 23
def istft_single_channel(stft_signal, size=1024, shift=256,
use_amplitude_for_biorthogonal_window=False,
disable_sythesis_window=False):
"""
Calculated the inverse short time Fourier transform to exactly reconstruct
the time signal.

Notes:
Be careful if you make modifications in the frequency domain (e.g.
beamforming) because the synthesis window is calculated according to the
unmodified! analysis window.

Args:
stft_signal: Single channel complex STFT signal
with dimensions frames times size/2+1.
size: Scalar FFT-size.
shift: Scalar FFT-shift. Typically shift is a fraction of size.
window: Window function handle.
window_length: Sometimes one desires to use a shorter window than
the fft size. In that case, the window is padded with zeros.
The default is to use the fft-size as a window size.

Returns:
Single channel complex STFT signal
Single channel time signal.
"""
assert stft_signal.shape[1] == size // 2 + 1, str(stft_signal.shape)

if window_length is None:
window = window(size)
else:
window = window(window_length)
window = np.pad(window, (0, size-window_length), mode='constant')
window = _biorthogonal_window_fastest(window, shift,
use_amplitude_for_biorthogonal_window)
if disable_sythesis_window:
window = np.ones_like(window)

time_signal = scipy.zeros(stft_signal.shape[0] * shift + size - shift)

for j, i in enumerate(range(0, len(time_signal) - size + shift, shift)):
time_signal[i:i + size] += window * np.real(irfft(stft_signal[j]))

time_signal = time_signal[size-shift:len(time_signal)-(size-shift)]

return time_signal 
Example 24
def roll_invariant_euclidean_distances(X, Y=None, squared=False):
"""
Considering the rows of X (and Y=X) as vectors, compute the
distance matrix between each pair of vectors.
The distance is the minimum of the euclidean distance over all rolls:

dist(x, y) = min_\tau(||x(t) - y(t - \tau)||^2)

Parameters
----------
X : array, shape (n_samples_1, n_features)

Y : array, shape (n_samples_2, n_features)

squared : boolean
Not used. Only for API compatibility.

Returns
-------
distances : array, shape (n_samples_1, n_samples_2)

"""
X = np.atleast_2d(X)
if Y is not None:
Y = np.atleast_2d(Y)
X, Y = check_pairwise_arrays(X, Y)
n_samples_1, n_features = X.shape
n_samples_2, n_features = Y.shape

X_norm = np.power(np.linalg.norm(X, axis=1), 2)
Y_norm = np.power(np.linalg.norm(Y, axis=1), 2)

# n_fft = next_fast_len(n_features + n_pads)
n_fft = n_features  # not fast but otherwise the distance is wrong
X_hat = rfft(X, n_fft, axis=1)
Y_hat = rfft(Y, n_fft, axis=1).conj()

# # broadcasting can have a huge memory cost
# XY_hat = X_hat[:, None, :] * Y_hat[None, :, :]
# XY = irfft(XY_hat, n_fft, axis=2).max(axis=2)
# distances = X_norm[:, None] + Y_norm[None, :] - 2 * XY

distances = np.zeros((n_samples_1, n_samples_2))
if n_samples_2 > 1:
print('RIED on %s samples, this might be slow' % (distances.shape, ))
for ii in range(n_samples_1):
for jj in range(n_samples_2):
XY = irfft(X_hat[ii] * Y_hat[jj], n_fft).max()
distances[ii, jj] = X_norm[ii] + Y_norm[jj] - 2 * XY

distances += 1e-12

return distances 
Example 25
def falpha(length=8192, alpha=1.0, fl=None, fu=None, mean=0.0, var=1.0):
"""Generate (1/f)^alpha noise by inverting the power spectrum.

Generates (1/f)^alpha noise by inverting the power spectrum.
Follows the algorithm described by Voss (1988) to generate
fractional Brownian motion.

Parameters
----------
length : int, optional (default = 8192)
Length of the time series to be generated.
alpha : float, optional (default = 1.0)
Exponent in (1/f)^alpha.  Pink noise will be generated by
default.
fl : float, optional (default = None)
Lower cutoff frequency.
fu : float, optional (default = None)
Upper cutoff frequency.
mean : float, optional (default = 0.0)
Mean of the generated noise.
var : float, optional (default = 1.0)
Variance of the generated noise.

Returns
-------
x : array
Array containing the time series.

Notes
-----
As discrete Fourier transforms assume that the input data is
periodic, the resultant series x_{i} (= x_{i + N}) is also periodic.
To avoid this periodicity, it is recommended to always generate
a longer series (two or three times longer) and trim it to the
desired length.
"""
freqs = fft.rfftfreq(length)
power = freqs[1:] ** -alpha
power = np.insert(power, 0, 0)  # P(0) = 0

if fl:
power[freqs < fl] = 0

if fu:
power[freqs > fu] = 0

# Randomize complex phases.
phase = 2 * np.pi * np.random.random(len(freqs))
y = np.sqrt(power) * np.exp(1j * phase)

# The last component (corresponding to the Nyquist frequency) of an
# RFFT with even number of points is always real.  (We don't have to
# make the mean real as P(0) = 0.)
if length % 2 == 0:
y[-1] = np.abs(y[-1] * np.sqrt(2))

x = fft.irfft(y, n=length)

# Rescale to proper variance and mean.
x = np.sqrt(var) * x / np.std(x)
return mean + x - np.mean(x)