# Python numpy.fft.rfft() Examples

The following are 23 code examples for showing how to use numpy.fft.rfft(). 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 bench_random(self):
from numpy.fft import rfft as numpy_rfft
print()
print('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)
print('|%8.2f' % measure('rfft(x)',repeat), end=' ')
sys.stdout.flush()

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

print(' (secs for %s calls)' % (repeat))
sys.stdout.flush() 
Example 4
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 5
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 6
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 7
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 8
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 9
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 10
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 11
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 12
def PowerSpectrum(X, Fs):
"""
Calculate the power spectrum of real-valued signal X measured with a
sampling frequency Fs. Result is in the same units as Fs.
"""

X = X - X.mean()
F = rfft(X)[:-1]
freq = np.arange(0, 0.5*Fs, Fs*1.0/len(X))
pw = np.abs(F)**2
return freq, pw/max(pw) 
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 analysis(self, x):
"""
Perform frequency analysis of a real input using DFT.

Parameters
----------
x : numpy array
Real signal in time domain.

Returns
-------
numpy array
DFT of input.
"""

# check for valid input
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]))

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

# apply DFT
if self.transform == 'fftw':
self.a[:,] = x
self.X[:,] = self._forward()
elif self.transform == 'mkl':
self.X[:,] = mkl_fft.rfft_numpy(x, self.nfft, axis=self.axis)
else:
self.X[:,] = rfft(x, self.nfft, axis=self.axis)

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 processFrames(self, windowedSamples):
self.complexMixtureSpectrogram[:] = rfft(windowedSamples * self.windowFunction, axis=1).astype(np.complex64)
self.spectrogram.set_value(self.complexMixtureSpectrogram)
#self.spectrogram.set_value( rfft(windowedSamples * self.windowFunction, axis=1).astype(np.complex64) )

realGCC = self.getComplexGCC()[0].real
if self.separationEnabled:

else:
outputSpectrogram = self.complexMixtureSpectrogram.copy()

if self.inputSpectrogramHistory:
self.inputSpectrogramHistory.set( -np.mean(np.abs(self.complexMixtureSpectrogram), axis=0) ** (1/3.0) )
if self.gccPHATHistory:
self.gccPHATHistory.set( np.nanmean(realGCC, axis=0).T )
if self.tdoaHistory:
if self.localizationEnabled:
gccPHATHistory = self.gccPHATHistory.getUnraveledArray()
tdoaIndex = np.argmax( np.nanmean(gccPHATHistory[:, -self.localizationWindowSize:], axis=-1) )
#tdoaIndex = (self.targetTDOAIndex.get_value() + 1) % self.numTDOAs
#tdoaIndex = np.random.randint(0, self.numTDOAs+1)
self.targetTDOAIndex.set_value(tdoaIndex)
self.tdoaHistory.set( np.array( [[self.targetTDOAIndex.get_value()]] ) )
if self.outputSpectrogramHistory:
self.outputSpectrogramHistory.set( -np.nanmean(np.abs(outputSpectrogram), axis=0) ** (1/3.0) )

return np.fft.irfft(outputSpectrogram, axis=1) * self.synthesisWindowFunction 
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 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