我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用scipy.fftpack.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 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 frft(x, alpha): """ Implementation of the Fractional Fourier Transform (FRFT) :param x: array of data to be transformed :param alpha: parameter of FRFT :return: FRFT(x) """ k = np.arange(x.size) y = np.hstack([ x * np.exp(-np.pi * 1j * k**2 * alpha), np.zeros(x.size, dtype=np.complex) ]) z = np.hstack([ np.exp(np.pi * 1j * k**2 * alpha), np.exp(np.pi * 1j * (k - x.size)**2 * alpha) ]) G = fftpack.ifft( fftpack.fft(y, overwrite_x=True) * fftpack.fft(z, overwrite_x=True), overwrite_x=True ) return np.exp(-np.pi * 1j * k**2 * alpha) * G[:x.size] # generate the desired momentum grid
def spectral_whitening(data, sr=None, smooth=None, filter=None, waterlevel=1e-8): """ Apply spectral whitening to data Data is divided by its smoothed (Default: None) amplitude spectrum. :param data: numpy array with data to manipulate :param sr: sampling rate (only needed for smoothing) :param smooth: length of smoothing window in Hz (default None -> no smoothing) :param filter: filter spectrum with bandpass after whitening (tuple with min and max frequency) :param waterlevel: waterlevel relative to mean of spectrum :return: whitened data """ mask = np.ma.getmask(data) N = len(data) nfft = next_fast_len(N) spec = fft(data, nfft) spec_ampl = np.sqrt(np.abs(np.multiply(spec, np.conjugate(spec)))) if smooth: smooth = int(smooth * N / sr) spec_ampl = ifftshift(smooth_func(fftshift(spec_ampl), smooth)) # save guard against division by 0 wl = waterlevel * np.mean(spec_ampl) spec_ampl[spec_ampl < wl] = wl spec /= spec_ampl if filter is not None: spec *= _filter_resp(*filter, sr=sr, N=len(spec), whole=True)[1] ret = np.real(ifft(spec, nfft)[:N]) return _fill_array(ret, mask=mask, fill_value=0.)
def fft(fEM, time, freq, ftarg): """Fourier Transform using the Fast Fourier Transform. The function is called from one of the modelling routines in :mod:`model`. Consult these modelling routines for a description of the input and output parameters. Returns ------- tEM : array Returns time-domain EM response of `fEM` for given `time`. conv : bool Only relevant for QWE/QUAD. """ # Get ftarg values dfreq, nfreq, ntot, pts_per_dec = ftarg # If pts_per_dec, we have first to interpolate fEM to required freqs if pts_per_dec: sfEMr = iuSpline(np.log10(freq), fEM.real) sfEMi = iuSpline(np.log10(freq), fEM.imag) freq = np.arange(1, nfreq+1)*dfreq fEM = sfEMr(np.log10(freq)) + 1j*sfEMi(np.log10(freq)) # Pad the frequency result fEM = np.pad(fEM, (0, ntot-nfreq), 'linear_ramp') # Carry out FFT ifftEM = fftpack.ifft(np.r_[fEM[1:], 0, fEM[::-1].conj()]).real stEM = 2*ntot*fftpack.fftshift(ifftEM*dfreq, 0) # Interpolate in time domain dt = 1/(2*ntot*dfreq) ifEM = iuSpline(np.linspace(-ntot, ntot-1, 2*ntot)*dt, stEM) tEM = ifEM(time)/2*np.pi # (Multiplication of 2/pi in model.tem) # Return the electromagnetic time domain field # (Second argument is only for QWE) return tEM, True # 3. Utilities
def iAugFFT(Xaug,axis=0): F=Xaug.shape[axis]/2 X=np.take(Xaug,np.arange(0,F),axis=axis)+np.complex64(1j)*np.take(Xaug,np.arange(F,2*F),axis=axis) X=np.concatenate((X.conj(), np.take(X,np.arange(F-2,0,-1),axis=axis)), axis=axis) xr=fft.ifft(X,axis=axis).real return xr
def estimate(self, nbcorr=np.nan, numpsd=-1): fft_length, _ = self.check_params() if np.isnan((nbcorr)): nbcorr = self.ordar # -------- estimate correlation from psd full_psd = self.psd[numpsd] full_psd = np.c_[full_psd, np.conjugate(full_psd[:, :0:-1])] correl = fftpack.ifft(full_psd[0], fft_length, 0).real # -------- estimate AR part col1 = correl[self.ordma:self.ordma + nbcorr] row1 = correl[np.abs( np.arange(self.ordma, self.ordma - self.ordar, -1))] R = linalg.toeplitz(col1, row1) r = -correl[self.ordma + 1:self.ordma + nbcorr + 1] AR = linalg.solve(R, r) self.AR_ = AR # -------- estimate correlation of MA part # -------- estimate MA part if self.ordma == 0: sigma2 = correl[0] + np.dot(AR, correl[1:self.ordar + 1]) self.MA = np.ones(1) * np.sqrt(sigma2) else: raise NotImplementedError( 'arma: estimation of the MA part not yet implemented')
def fft_multiply_repeated(h_fft, x, cuda_dict=dict(use_cuda=False)): """Do FFT multiplication by a filter function (possibly using CUDA) Parameters ---------- h_fft : 1-d array or gpuarray The filtering array to apply. x : 1-d array The array to filter. cuda_dict : dict Dictionary constructed using setup_cuda_multiply_repeated(). Returns ------- x : 1-d array Filtered version of x. """ if not cuda_dict['use_cuda']: # do the fourier-domain operations x = np.real(ifft(h_fft * fft(x), overwrite_x=True)).ravel() else: cudafft = _get_cudafft() # do the fourier-domain operations, results in second param cuda_dict['x'].set(x.astype(np.float64)) cudafft.fft(cuda_dict['x'], cuda_dict['x_fft'], cuda_dict['fft_plan']) _multiply_inplace_c128(h_fft, cuda_dict['x_fft']) # If we wanted to do it locally instead of using our own kernel: # cuda_seg_fft.set(cuda_seg_fft.get() * h_fft) cudafft.ifft(cuda_dict['x_fft'], cuda_dict['x'], cuda_dict['ifft_plan'], False) x = np.array(cuda_dict['x'].get(), dtype=x.dtype, subok=True, copy=False) return x ############################################################################### # FFT Resampling
def _st(x, start_f, windows): """Implementation based on Ali Moukadem Matlab code (only used in tests)""" n_samp = x.shape[-1] ST = np.empty(x.shape[:-1] + (len(windows), n_samp), dtype=np.complex) # do the work Fx = fftpack.fft(x) XF = np.concatenate([Fx, Fx], axis=-1) for i_f, window in enumerate(windows): f = start_f + i_f ST[..., i_f, :] = fftpack.ifft(XF[..., f:f + n_samp] * window) return ST
def shift_waveform(tr, dtshift): """ tr_shift = shift_waveform(tr, dtshift): Shift data in trace tr by dtshift seconds backwards. Parameters ---------- tr : obspy.Trace Trace that contains the data to shift dtshift : float Time shift in seconds Returns ------- tr_shift : obspy.Trace Copy of tr, but with data shifted dtshift seconds backwards. """ data_pad = np.r_[tr.data, np.zeros_like(tr.data)] freq = fft.fftfreq(len(data_pad), tr.stats.delta) shiftvec = np.exp(- 2 * np.pi * complex(0., 1.) * freq * dtshift) data_fd = shiftvec * fft.fft(data_pad * signal.tukey(len(data_pad), alpha=0.2)) tr.data = np.real(fft.ifft(data_fd))[0:tr.stats.npts]
def phase_randomize(D, random_state=0): """Randomly shift signal phases For each timecourse (from each voxel and each subject), computes its DFT and then randomly shifts the phase of each frequency before inverting back into the time domain. This yields timecourses with the same power spectrum (and thus the same autocorrelation) as the original timecourses, but will remove any meaningful temporal relationships between the timecourses. This procedure is described in: Simony E, Honey CJ, Chen J, Lositsky O, Yeshurun Y, Wiesel A, Hasson U (2016) Dynamic reconfiguration of the default mode network during narrative comprehension. Nat Commun 7. Parameters ---------- D : voxel by time by subject ndarray fMRI data to be phase randomized random_state : RandomState or an int seed (0 by default) A random number generator instance to define the state of the random permutations generator. Returns ---------- ndarray of same shape as D phase randomized timecourses """ random_state = check_random_state(random_state) F = fft(D, axis=1) if D.shape[1] % 2 == 0: pos_freq = np.arange(1, D.shape[1] // 2) neg_freq = np.arange(D.shape[1] - 1, D.shape[1] // 2, -1) else: pos_freq = np.arange(1, (D.shape[1] - 1) // 2 + 1) neg_freq = np.arange(D.shape[1] - 1, (D.shape[1] - 1) // 2, -1) shift = random_state.rand(D.shape[0], len(pos_freq), D.shape[2]) * 2 * math.pi # Shift pos and neg frequencies symmetrically, to keep signal real F[:, pos_freq, :] *= np.exp(1j * shift) F[:, neg_freq, :] *= np.exp(-1j * shift) return np.real(ifft(F, axis=1))
def nsgitf_real(c, c_dc, c_nyq, multiscale, shift): """ Nonstationary Inverse Gabor Transform on real valued signal References ---------- Velasco G. A., Holighaus N., Dorfler M., Grill T. Constructing an invertible constant-Q transform with nonstationary Gabor frames, Proceedings of the 14th International Conference on Digital Audio Effects (DAFx 11), Paris, France, 2011 Holighaus N., Dorfler M., Velasco G. A. and Grill T. A framework for invertible, real-time constant-Q transforms, submitted. Original matlab code copyright follows: AUTHOR(s) : Monika Dorfler, Gino Angelo Velasco, Nicki Holighaus, 2010-2011 COPYRIGHT : (c) NUHAG, Dept.Math., University of Vienna, AUSTRIA http://nuhag.eu/ Permission is granted to modify and re-distribute this code in any manner as long as this notice is preserved. All standard disclaimers apply. """ c_l = [] c_l.append(c_dc) c_l.extend([ci for ci in c]) c_l.append(c_nyq) posit = np.cumsum(shift) seq_len = posit[-1] posit -= shift[0] out = np.zeros((seq_len,)).astype(c_l[1].dtype) for ii in range(len(c_l)): filt_len = len(multiscale[ii]) win_range = posit[ii] + np.arange(-np.floor(filt_len / 2), np.ceil(filt_len / 2)) win_range = (win_range % seq_len).astype(np.int) temp = np.fft.fft(c_l[ii]) * len(c_l[ii]) fs_new_bins = len(c_l[ii]) fk_bins = posit[ii] displace = int(fk_bins - np.floor(fk_bins / fs_new_bins) * fs_new_bins) temp = np.roll(temp, -displace) l = np.arange(len(temp) - np.floor(filt_len / 2), len(temp)) r = np.arange(np.ceil(filt_len / 2)) temp_idx = (np.concatenate((l, r)) % len(temp)).astype(np.int) temp = temp[temp_idx] lf = np.arange(filt_len - np.floor(filt_len / 2), filt_len) rf = np.arange(np.ceil(filt_len / 2)) filt_idx = np.concatenate((lf, rf)).astype(np.int) m = multiscale[ii][filt_idx] out[win_range] = out[win_range] + m * temp nyq_bin = np.floor(seq_len / 2) + 1 out_idx = np.arange( nyq_bin - np.abs(1 - seq_len % 2) - 1, 0, -1).astype(np.int) out[nyq_bin:] = np.conj(out[out_idx]) t_out = np.real(np.fft.ifft(out)).astype(np.float64) return t_out
def invert_spectrogram(X_s, step, calculate_offset=True, set_zero_phase=True): """ Under MSR-LA License Based on MATLAB implementation from Spectrogram Inversion Toolbox References ---------- D. Griffin and J. Lim. Signal estimation from modified short-time Fourier transform. IEEE Trans. Acoust. Speech Signal Process., 32(2):236-243, 1984. Malcolm Slaney, Daniel Naar and Richard F. Lyon. Auditory Model Inversion for Sound Separation. Proc. IEEE-ICASSP, Adelaide, 1994, II.77-80. Xinglei Zhu, G. Beauregard, L. Wyse. Real-Time Signal Estimation from Modified Short-Time Fourier Transform Magnitude Spectra. IEEE Transactions on Audio Speech and Language Processing, 08/2007. """ size = int(X_s.shape[1] // 2) wave = np.zeros((X_s.shape[0] * step + size)) # Getting overflow warnings with 32 bit... wave = wave.astype('float64') total_windowing_sum = np.zeros((X_s.shape[0] * step + size)) win = 0.54 - .46 * np.cos(2 * np.pi * np.arange(size) / (size - 1)) est_start = int(size // 2) - 1 est_end = est_start + size for i in range(X_s.shape[0]): wave_start = int(step * i) wave_end = wave_start + size if set_zero_phase: spectral_slice = X_s[i].real + 0j else: # already complex spectral_slice = X_s[i] # Don't need fftshift due to different impl. wave_est = np.real(np.fft.ifft(spectral_slice))[::-1] if calculate_offset and i > 0: offset_size = size - step if offset_size <= 0: print("WARNING: Large step size >50\% detected! " "This code works best with high overlap - try " "with 75% or greater") offset_size = step offset = xcorr_offset(wave[wave_start:wave_start + offset_size], wave_est[est_start:est_start + offset_size]) else: offset = 0 wave[wave_start:wave_end] += win * wave_est[ est_start - offset:est_end - offset] total_windowing_sum[wave_start:wave_end] += win wave = np.real(wave) / (total_windowing_sum + 1E-6) return wave
def test_ffts(self): t, tsc, y, err = data() yhat = np.empty(len(y)) yg = gpuarray.to_gpu(y.astype(np.complex128)) yghat = gpuarray.to_gpu(yhat.astype(np.complex128)) plan = cufft.Plan(len(y), np.complex128, np.complex128) cufft.ifft(yg, yghat, plan) yhat = fftpack.ifft(y) * len(y) tols = dict(rtol=nfft_rtol, atol=nfft_atol) assert_allclose(yhat, yghat.get(), **tols)
def fftconvolve(in1, in2, mode="full", axis=None): """ Convolve two N-dimensional arrays using FFT. See convolve. This is a fix of scipy.signal.fftconvolve, adding an axis argument and importing locally the stuff only needed for this function """ s1 = np.array(in1.shape) s2 = np.array(in2.shape) complex_result = (np.issubdtype(in1.dtype, np.complex) or np.issubdtype(in2.dtype, np.complex)) if axis is None: size = s1 + s2 - 1 fslice = tuple([slice(0, int(sz)) for sz in size]) else: equal_shapes = s1 == s2 # allow equal_shapes[axis] to be False equal_shapes[axis] = True assert equal_shapes.all(), 'Shape mismatch on non-convolving axes' size = s1[axis] + s2[axis] - 1 fslice = [slice(l) for l in s1] fslice[axis] = slice(0, int(size)) fslice = tuple(fslice) # Always use 2**n-sized FFT fsize = 2 ** int(np.ceil(np.log2(size))) if axis is None: IN1 = fftpack.fftn(in1, fsize) IN1 *= fftpack.fftn(in2, fsize) ret = fftpack.ifftn(IN1)[fslice].copy() else: IN1 = fftpack.fft(in1, fsize, axis=axis) IN1 *= fftpack.fft(in2, fsize, axis=axis) ret = fftpack.ifft(IN1, axis=axis)[fslice].copy() del IN1 if not complex_result: ret = ret.real if mode == "full": return ret elif mode == "same": if np.product(s1, axis=0) > np.product(s2, axis=0): osize = s1 else: osize = s2 return signaltools._centered(ret, osize) elif mode == "valid": return signaltools._centered(ret, abs(s2 - s1) + 1)
