def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2 + 1] = X
X_pad[:, fftsize // 2 + 1:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
python类ifft()的实例源码
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2 + 1] = X
X_pad[:, fftsize // 2 + 1:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def carr_madan_fraction_fft_call_pricer(N, d_u, d_k, alpha, r, t, S0, q, chf_ln_st):
rou = (d_u * d_k) / (2 * np.pi)
beta = np.log(S0) - d_k * N / 2
u_arr = np.arange(N) * d_u
k_arr = beta + np.arange(N) * d_k
delta_arr = np.zeros(N)
delta_arr[0] = 1
w_arr = d_u / 3 * (3 + (-1) ** (np.arange(N) + 1) - delta_arr)
call_chf = (np.exp(-r * t) / ((alpha + 1j * u_arr) * (alpha + 1j * u_arr + 1))) * chf_ln_st(
u_arr - (alpha + 1) * 1j,
t, r, q=q, S0=S0)
x_arr = np.exp(-1j * beta * u_arr) * call_chf * w_arr
y_arr = np.zeros(2 * N) * 0j
y_arr[:N] = np.exp(-1j * np.pi * rou * np.arange(N) ** 2) * x_arr
z_arr = np.zeros(2 * N) * 0j
z_arr[:N] = np.exp(1j * np.pi * rou * np.arange(N) ** 2)
z_arr[N:] = np.exp(1j * np.pi * rou * np.arange(N - 1, -1, -1) ** 2)
ffty = (fft(y_arr))
fftz = (fft(z_arr))
fftx = ffty * fftz
fftpsi = ifft(fftx)
fft_prices = np.exp(-1j * np.pi * (np.arange(N) ** 2) * rou) * fftpsi[:N]
call_prices = (np.exp(-alpha * k_arr) / np.pi) * fft_prices.real
return np.exp(k_arr), call_prices
def _lpc(x,order):
"""Linear predictor coefficients. Supports 1D and 2D arrays only through
iteration."""
x = np.asarray(x)
if x.ndim == 1:
m = x.size
X = fft(x,n=nextpow2(m)) # TODO nextpower of 2
R = np.real(ifft(np.abs(X)**2)) # Auto-correlation matrix
R = R/m
a = _levinson(r=R, order=order)[0]
a = np.hstack((1., a))
elif x.ndim == 2:
m,n = x.shape
X = fft(x,n=nextpow2(m), axis=0)
R = np.real(ifft(np.abs(X)**2, axis=0)) # Auto-correlation matrix
R = R/m
a = np.ones((order+1,n))
for col in range(n):
a[1:,col] = _levinson(r=R[:,col], order=order)[0]
else:
raise ValueError('Supported for 1-D or 2-D arrays only.')
return a
def compare_fft(date, x, file_name='comparefft.jpg'):
from scipy.fftpack import fft, ifft
ax = []
l = len(x)
t = 50
x = exp_moving_average(x, 10)
labels = []
ax.append(x[t:l-t])
labels.append("??????? ???")
ax.append((fft(x)[t:l-t]))
labels.append("???’? ????????????")
# print(l, len(np.fft.ifft(np.fft.fft(x)[5:l-5])[5:l-5]))
print(len(x))
# print(len(ifft(np.log(np.abs(fft(x)[5:l-5])))[5:l-5]))
# ax.append((ifft(np.log(np.abs(fft(x))))[t:l-t])[:l/2 - t])
# at.append((date[t:l-t])[::2])
print("---")
showPlotLabelsCompare(ax, date[t:l-t], labels, file_name)
def compare_cepstrum(date, x, file_name='comparecep.jpg'):
# date, x = hist_data.get_historical_quotes(start_date=datetime.datetime(2000, 2, 2),
# end_date=datetime.datetime(2001, 2, 2),
# csv_path='../../data/usdgbp1990.csv')
from scipy.fftpack import fft, ifft
ax = []
at = []
l = len(x)
t = 50
x = exp_moving_average(x, 10)
ax.append(x[t:l-t])
at.append(date[t:l-t])
# print(l, len(np.fft.ifft(np.fft.fft(x)[5:l-5])[5:l-5]))
print(len(x))
# print(len(ifft(np.log(np.abs(fft(x)[5:l-5])))[5:l-5]))
ax.append((ifft(np.log(np.abs(fft(x))))[t:l-t])[:l/2 - t])
at.append((date[t:l-t])[::2])
print("---")
showPlotCompareSeparate(ax, at, file_name)
def _st_power_itc(x, start_f, compute_itc, zero_pad, decim, W):
"""Aux function"""
n_samp = x.shape[-1]
n_out = (n_samp - zero_pad)
n_out = n_out // decim + bool(n_out % decim)
psd = np.empty((len(W), n_out))
itc = np.empty_like(psd) if compute_itc else None
X = fftpack.fft(x)
XX = np.concatenate([X, X], axis=-1)
for i_f, window in enumerate(W):
f = start_f + i_f
ST = fftpack.ifft(XX[:, f:f + n_samp] * window)
TFR = ST[:, :-zero_pad:decim]
TFR_abs = np.abs(TFR)
if compute_itc:
TFR /= TFR_abs
itc[i_f] = np.abs(np.mean(TFR, axis=0))
TFR_abs *= TFR_abs
psd[i_f] = np.mean(TFR_abs, axis=0)
return psd, itc
def test_pdos_1d():
pad=lambda x: pad_zeros(x, nadd=len(x)-1)
n=500; w=welch(n)
# 1 second signal
t=np.linspace(0,1,n); dt=t[1]-t[0]
# sum of sin()s with random freq and phase shift, 10 frequencies from
# f=0...100 Hz
v=np.array([np.sin(2*np.pi*f*t + rand()*2*np.pi) for f in rand(10)*100]).sum(0)
f=np.fft.fftfreq(2*n-1, dt)[:n]
c1=mirror(ifft(abs(fft(pad(v)))**2.0)[:n].real)
c2=correlate(v,v,'full')
c3=mirror(acorr(v,norm=False))
assert np.allclose(c1, c2)
assert np.allclose(c1, c3)
p1=(abs(fft(pad(v)))**2.0)[:n]
p2=(abs(fft(mirror(acorr(v,norm=False)))))[:n]
assert np.allclose(p1, p2)
p1=(abs(fft(pad(v*w)))**2.0)[:n]
p2=(abs(fft(mirror(acorr(v*w,norm=False)))))[:n]
assert np.allclose(p1, p2)
def polyval(self, chebcoeff):
"""
Compute the interpolation values at Chebyshev points.
chebcoeff: Chebyshev coefficients
"""
N = len(chebcoeff)
if N == 1:
return chebcoeff
data = even_data(chebcoeff)/2
data[0] *= 2
data[N-1] *= 2
fftdata = 2*(N-1)*fftpack.ifft(data, axis=0)
complex_values = fftdata[:N]
# convert to real if input was real
if np.isrealobj(chebcoeff):
values = np.real(complex_values)
else:
values = complex_values
return values
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def istft(X, fftsize=128, step="half", wsola=False, mean_normalize=True,
real=False, compute_onesided=True):
"""
Compute ISTFT for STFT transformed X
"""
if real:
local_ifft = fftpack.irfft
X_pad = np.zeros((X.shape[0], X.shape[1] + 1)) + 0j
X_pad[:, :-1] = X
X = X_pad
else:
local_ifft = fftpack.ifft
if compute_onesided:
X_pad = np.zeros((X.shape[0], 2 * X.shape[1])) + 0j
X_pad[:, :fftsize // 2] = X
X_pad[:, fftsize // 2:] = 0
X = X_pad
X = local_ifft(X).astype("float64")
if step == "half":
X = invert_halfoverlap(X)
else:
X = overlap_add(X, step, wsola=wsola)
if mean_normalize:
X -= np.mean(X)
return X
def iDFT(magX, phsX, wsz):
""" Discrete Fourier Transformation(Synthesis) of a given spectral analysis
via an inverse FFT implementation from scipy.
Args:
magX : (2D ndarray) Magnitude Spectrum
phsX : (2D ndarray) Phase Spectrum
wsz : (int) Synthesis window size
Returns:
y : (array) Real time domain output signal
"""
# Get FFT Size
hlfN = magX.size
N = (hlfN-1)*2
# Half of window size parameters
hw1 = int(math.floor((wsz+1)/2))
hw2 = int(math.floor(wsz/2))
# Initialise output spectrum with zeros
Y = np.zeros(N, dtype=complex)
# Initialise output array with zeros
y = np.zeros(wsz)
# Compute complex spectrum(both sides) in two steps
Y[0:hlfN] = magX * np.exp(1j*phsX)
Y[hlfN:] = magX[-2:0:-1] * np.exp(-1j*phsX[-2:0:-1])
# Perform the iDFT
fftbuffer = np.real(ifft(Y))
# Roll-back the zero-phase windowing technique
y[:hw2] = fftbuffer[-hw2:]
y[hw2:] = fftbuffer[:hw1]
return y
def iDFT(magX, phsX, wsz):
""" Discrete Fourier Transformation(Synthesis) of a given spectral analysis
via an inverse FFT implementation from scipy.
Args:
magX : (2D ndarray) Magnitude Spectrum
phsX : (2D ndarray) Phase Spectrum
wsz : (int) Synthesis window size
Returns:
y : (array) Real time domain output signal
"""
# Get FFT Size
hlfN = magX.size;
N = (hlfN-1)*2
# Half of window size parameters
hw1 = int(math.floor((wsz+1)/2))
hw2 = int(math.floor(wsz/2))
# Initialise synthesis buffer with zeros
fftbuffer = np.zeros(N)
# Initialise output spectrum with zeros
Y = np.zeros(N, dtype = complex)
# Initialise output array with zeros
y = np.zeros(wsz)
# Compute complex spectrum(both sides) in two steps
Y[0:hlfN] = magX * np.exp(1j*phsX)
Y[hlfN:] = magX[-2:0:-1] * np.exp(-1j*phsX[-2:0:-1])
# Perform the iDFT
fftbuffer = np.real(ifft(Y))
# Roll-back the zero-phase windowing technique
y[:hw2] = fftbuffer[-hw2:]
y[hw2:] = fftbuffer[:hw1]
return y
def _shift_fft(array, shift_value):
Ndim = array.ndim
dims = array.shape
dtype = array.dtype.kind
if (dtype != 'f'):
raise ValueError('Array must be float')
shifted = array
if (Ndim == 1):
Nx = dims[0]
x_ramp = np.arange(Nx, dtype=array.dtype) - Nx//2
tilt = (2*np.pi/Nx) * (shift_value[0]*x_ramp)
cplx_tilt = np.cos(tilt) + 1j*np.sin(tilt)
cplx_tilt = fft.fftshift(cplx_tilt)
narray = fft.fft(fft.ifft(array) * cplx_tilt)
shifted = narray.real
elif (Ndim == 2):
Nx = dims[0]
Ny = dims[1]
x_ramp = np.outer(np.full(Nx, 1.), np.arange(Ny, dtype=array.dtype)) - Nx//2
y_ramp = np.outer(np.arange(Nx, dtype=array.dtype), np.full(Ny, 1.)) - Ny//2
tilt = (2*np.pi/Nx) * (shift_value[0]*x_ramp+shift_value[1]*y_ramp)
cplx_tilt = np.cos(tilt) + 1j*np.sin(tilt)
cplx_tilt = fft.fftshift(cplx_tilt)
narray = fft.fft2(fft.ifft2(array) * cplx_tilt)
shifted = narray.real
else:
raise ValueError('This function can shift only 1D or 2D arrays')
return shifted
def harvest_get_f0_candidate_from_raw_event(boundary_f0,
fs, y_spectrum, y_length, temporal_positions, f0_floor,
f0_ceil):
filter_length_half = int(np.round(fs / boundary_f0 * 2))
band_pass_filter_base = harvest_nuttall(filter_length_half * 2 + 1)
shifter = np.cos(2 * np.pi * boundary_f0 * np.arange(-filter_length_half, filter_length_half + 1) / float(fs))
band_pass_filter = band_pass_filter_base * shifter
index_bias = filter_length_half
# possible numerical issues if 32 bit
spectrum_low_pass_filter = np.fft.fft(band_pass_filter.astype("float64"), len(y_spectrum))
filtered_signal = np.real(np.fft.ifft(spectrum_low_pass_filter * y_spectrum))
index_bias = filter_length_half + 1
filtered_signal = filtered_signal[index_bias + np.arange(y_length).astype("int32")]
negative_zero_cross = harvest_zero_crossing_engine(filtered_signal, fs)
positive_zero_cross = harvest_zero_crossing_engine(-filtered_signal, fs)
d_filtered_signal = filtered_signal[1:] - filtered_signal[:-1]
peak = harvest_zero_crossing_engine(d_filtered_signal, fs)
dip = harvest_zero_crossing_engine(-d_filtered_signal, fs)
f0_candidate = harvest_get_f0_candidate_contour(negative_zero_cross,
positive_zero_cross, peak, dip, temporal_positions)
f0_candidate[f0_candidate > (boundary_f0 * 1.1)] = 0.
f0_candidate[f0_candidate < (boundary_f0 * .9)] = 0.
f0_candidate[f0_candidate > f0_ceil] = 0.
f0_candidate[f0_candidate < f0_floor] = 0.
return f0_candidate
def cheaptrick_smoothing_with_recovery(smoothed_spectrum, f0, fs, fft_size, q1):
quefrency_axis = np.arange(fft_size) / float(fs)
# 0 is NaN
smoothing_lifter = np.sin(np.pi * f0 * quefrency_axis) / (np.pi * f0 * quefrency_axis)
p = smoothing_lifter[1:int(fft_size / 2)][::-1].copy()
smoothing_lifter[int(fft_size / 2) + 1:] = p
smoothing_lifter[0] = 1.
compensation_lifter = (1 - 2. * q1) + 2. * q1 * np.cos(2 * np.pi * quefrency_axis * f0)
p = compensation_lifter[1:int(fft_size / 2)][::-1].copy()
compensation_lifter[int(fft_size / 2) + 1:] = p
tandem_cepstrum = np.fft.fft(np.log(smoothed_spectrum))
tmp_spectral_envelope = np.exp(np.real(np.fft.ifft(tandem_cepstrum * smoothing_lifter * compensation_lifter)))
spectral_envelope = tmp_spectral_envelope[:int(fft_size / 2) + 1]
return spectral_envelope
def harvest_get_f0_candidate_from_raw_event(boundary_f0,
fs, y_spectrum, y_length, temporal_positions, f0_floor,
f0_ceil):
filter_length_half = int(np.round(fs / boundary_f0 * 2))
band_pass_filter_base = harvest_nuttall(filter_length_half * 2 + 1)
shifter = np.cos(2 * np.pi * boundary_f0 * np.arange(-filter_length_half, filter_length_half + 1) / float(fs))
band_pass_filter = band_pass_filter_base * shifter
index_bias = filter_length_half
# possible numerical issues if 32 bit
spectrum_low_pass_filter = np.fft.fft(band_pass_filter.astype("float64"), len(y_spectrum))
filtered_signal = np.real(np.fft.ifft(spectrum_low_pass_filter * y_spectrum))
index_bias = filter_length_half + 1
filtered_signal = filtered_signal[index_bias + np.arange(y_length).astype("int32")]
negative_zero_cross = harvest_zero_crossing_engine(filtered_signal, fs)
positive_zero_cross = harvest_zero_crossing_engine(-filtered_signal, fs)
d_filtered_signal = filtered_signal[1:] - filtered_signal[:-1]
peak = harvest_zero_crossing_engine(d_filtered_signal, fs)
dip = harvest_zero_crossing_engine(-d_filtered_signal, fs)
f0_candidate = harvest_get_f0_candidate_contour(negative_zero_cross,
positive_zero_cross, peak, dip, temporal_positions)
f0_candidate[f0_candidate > (boundary_f0 * 1.1)] = 0.
f0_candidate[f0_candidate < (boundary_f0 * .9)] = 0.
f0_candidate[f0_candidate > f0_ceil] = 0.
f0_candidate[f0_candidate < f0_floor] = 0.
return f0_candidate
def cheaptrick_smoothing_with_recovery(smoothed_spectrum, f0, fs, fft_size, q1):
quefrency_axis = np.arange(fft_size) / float(fs)
# 0 is NaN
smoothing_lifter = np.sin(np.pi * f0 * quefrency_axis) / (np.pi * f0 * quefrency_axis)
p = smoothing_lifter[1:int(fft_size / 2)][::-1].copy()
smoothing_lifter[int(fft_size / 2) + 1:] = p
smoothing_lifter[0] = 1.
compensation_lifter = (1 - 2. * q1) + 2. * q1 * np.cos(2 * np.pi * quefrency_axis * f0)
p = compensation_lifter[1:int(fft_size / 2)][::-1].copy()
compensation_lifter[int(fft_size / 2) + 1:] = p
tandem_cepstrum = np.fft.fft(np.log(smoothed_spectrum))
tmp_spectral_envelope = np.exp(np.real(np.fft.ifft(tandem_cepstrum * smoothing_lifter * compensation_lifter)))
spectral_envelope = tmp_spectral_envelope[:int(fft_size / 2) + 1]
return spectral_envelope
def iDFT(magX, phsX, wsz):
""" Discrete Fourier Transformation(Synthesis) of a given spectral analysis
via an inverse FFT implementation from scipy.
Args:
magX : (2D ndarray) Magnitude Spectrum
phsX : (2D ndarray) Phase Spectrum
wsz : (int) Synthesis window size
Returns:
y : (array) Real time domain output signal
"""
# Get FFT Size
hlfN = magX.size;
N = (hlfN-1)*2
# Half of window size parameters
hw1 = int(math.floor((wsz+1)/2))
hw2 = int(math.floor(wsz/2))
# Initialise synthesis buffer with zeros
fftbuffer = np.zeros(N)
# Initialise output spectrum with zeros
Y = np.zeros(N, dtype = complex)
# Initialise output array with zeros
y = np.zeros(wsz)
# Compute complex spectrum(both sides) in two steps
Y[0:hlfN] = magX * np.exp(1j*phsX)
Y[hlfN:] = magX[-2:0:-1] * np.exp(-1j*phsX[-2:0:-1])
# Perform the iDFT
fftbuffer = np.real(ifft(Y))
# Roll-back the zero-phase windowing technique
y[:hw2] = fftbuffer[-hw2:]
y[hw2:] = fftbuffer[:hw1]
return y
