The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
def crosscorr(tr1,tr2):
"""This function takes 2 obspy traces and performs a phase crosscorrelation
of the traces. First, the hilbert transform is taken to obtain the
analytic signal and hence the instantaneous phase. This is then passed to a
fortran script where phase correlation is performed after Schimmel et al., 2011.
This function relies on the shared object file phasecorr.so, which is the file
containing the fortran subroutine for phase correlation.
"""
import numpy as np
from scipy.signal import hilbert
from phasecorr import phasecorr
# Hilbert transform to obtain the analytic signal
htrans1 = hilbert(tr1)
htrans2 = hilbert(tr2)
# Perform phase autocorrelation with instantaneous phase
pcc = phasecorr(np.angle(htrans1),np.angle(htrans2),len(htrans1))
return pcc Example 2
def getAngleInDegrees(_fft):
_angle = np.angle(_fft, True);
_rows, _columns = _fft.shape;
for i in range(0, _rows):
for j in range(0, _columns):
if (_angle[i][j] < 0 and abs(_angle[i][j]) <= 180):
#print(_angle[i][j]);
#print("\n");
_angle[i][j] = 180 + _angle[i][j];
#print(_angle[i][j]);
#print("\n");
elif (_angle[i][j] < 0 and abs(_angle[i][j]) > 180):
#print(_angle[i][j]);
#print("\n");
_angle[i][j] = 360 + _angle[i][j];
#print(_angle[i][j]);
#print("\n");
return _angle; Example 3
def appres(F, H, sig, chg, taux, c, mu, eps, n):
Res = np.zeros_like(F)
Phase = np.zeros_like(F)
App_ImpZ= np.zeros_like(F, dtype='complex_')
for i in range(0, len(F)):
UD, EH, Z , K = Propagate(F[i], H, sig, chg, taux, c, mu, eps, n)
App_ImpZ[i] = EH[0, 1]/EH[1, 1]
Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
Phase[i] = np.angle(App_ImpZ[i], deg = True)
return Res, Phase
# Evaluate Up, Down components, E and H field, for a frequency range,
# a discretized depth range and a time range (use to calculate envelope) Example 4
def phase_plot(data, toffset, log_scale, title):
"""Plot the phase of the data in linear or log scale."""
print("phase")
phase = numpy.angle(data) / numpy.pi
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(phase)
axmx = numpy.max(phase)
#ax.axis([-axmx, axmx, -axmx, axmx])
ax.grid(True)
ax.set_xlabel('time')
ax.set_ylabel('phase')
ax.set_title(title)
return fig Example 5
def unknown_feature_extractor(x, sr, win_len, shift_len, barks, inner_win, inner_shift, win_type, method_version):
x_spectrum = stft_extractor(x, win_len, shift_len, win_type)
coef = get_fft_bark_mat(sr, win_len, barks, 20, sr//2)
bark_spect = np.matmul(coef, x_spectrum)
ams = np.zeros((barks, inner_win//2+1, (bark_spect.shape[1] - inner_win)//inner_shift))
for i in range(barks):
channel_stft = stft_extractor(bark_spect[i, :], inner_win, inner_shift, 'hanning')
if method_version == 'v1':
ams[i, :, :] = 20 * np.log(np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift]))
elif method_version == 'v2':
channel_amplitude = np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
channel_angle = np.angle(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
channel_angle = channel_angle - (np.floor(channel_angle / (2.*np.pi)) * (2.*np.pi))
ams[i, :, :] = np.power(channel_amplitude, 1./3.) * channel_angle
else:
ams[i, :, :] = np.abs(channel_stft)
return ams Example 6
def ams_extractor(x, sr, win_len, shift_len, barks, inner_win, inner_shift, win_type, method_version):
x_spectrum = stft_extractor(x, win_len, shift_len, win_type)
coef = get_fft_bark_mat(sr, win_len, barks, 20, sr//2)
bark_spect = np.matmul(coef, x_spectrum)
ams = np.zeros((barks, inner_win//2+1, (bark_spect.shape[1] - inner_win)//inner_shift))
for i in range(barks):
channel_stft = stft_extractor(bark_spect[i, :], inner_win, inner_shift, 'hanning')
if method_version == 'v1':
ams[i, :, :] = 20 * np.log(np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift]))
elif method_version == 'v2':
channel_amplitude = np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
channel_angle = np.angle(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
channel_angle = channel_angle - (np.floor(channel_angle / (2.*np.pi)) * (2.*np.pi))
ams[i, :, :] = np.power(channel_amplitude, 1./3.) * channel_angle
else:
ams[i, :, :] = np.abs(channel_stft)
return ams Example 7
def griffin_lim(mag, phase_angle, n_fft, hop, num_iters):
"""Iterative algorithm for phase retrival from a magnitude spectrogram.
Args:
mag: Magnitude spectrogram.
phase_angle: Initial condition for phase.
n_fft: Size of the FFT.
hop: Stride of FFT. Defaults to n_fft/2.
num_iters: Griffin-Lim iterations to perform.
Returns:
audio: 1-D array of float32 sound samples.
"""
fft_config = dict(n_fft=n_fft, win_length=n_fft, hop_length=hop, center=True)
ifft_config = dict(win_length=n_fft, hop_length=hop, center=True)
complex_specgram = inv_magphase(mag, phase_angle)
for i in range(num_iters):
audio = librosa.istft(complex_specgram, **ifft_config)
if i != num_iters - 1:
complex_specgram = librosa.stft(audio, **fft_config)
_, phase = librosa.magphase(complex_specgram)
phase_angle = np.angle(phase)
complex_specgram = inv_magphase(mag, phase_angle)
return audio Example 8
def __init__(self, qubit_names, quad="real"):
super(PulseCalibration, self).__init__()
self.qubit_names = qubit_names if isinstance(qubit_names, list) else [qubit_names]
self.qubit = [QubitFactory(qubit_name) for qubit_name in qubit_names] if isinstance(qubit_names, list) else QubitFactory(qubit_names)
self.filename = 'None'
self.exp = None
self.axis_descriptor = None
self.cw_mode = False
self.saved_settings = config.load_meas_file(config.meas_file)
self.settings = deepcopy(self.saved_settings) #make a copy for used during calibration
self.quad = quad
if quad == "real":
self.quad_fun = np.real
elif quad == "imag":
self.quad_fun = np.imag
elif quad == "amp":
self.quad_fun = np.abs
elif quad == "phase":
self.quad_fun = np.angle
else:
raise ValueError('Quadrature to calibrate must be one of ("real", "imag", "amp", "phase").')
self.plot = self.init_plot() Example 9
def displayResult(self):
fig = plt.figure(101)
plt.subplot(221)
plt.imshow(np.abs(self.reconstruction),origin='lower')
plt.draw()
plt.title('Reconstruction Magnitude')
plt.subplot(222)
plt.imshow(np.angle(self.reconstruction),origin='lower')
plt.draw()
plt.title('Reconstruction Phase')
plt.subplot(223)
plt.imshow(np.abs(self.aperture),origin='lower')
plt.title("Aperture Magnitude")
plt.draw()
plt.subplot(224)
plt.imshow(np.angle(self.aperture),origin='lower')
plt.title("Aperture Phase")
plt.draw()
fig.canvas.draw()
#fig.tight_layout()
# display.display(fig)
# display.clear_output(wait=True)
# time.sleep(.00001) Example 10
def displayResult2(self):
fig = plt.figure()
ax = fig.add_subplot(2,2,1 )
ax.imshow(np.abs(self.reconstruction))
ax.set_title('Reconstruction Magnitude')
ax = fig.add_subplot(2,2,2 )
ax.imshow(np.angle(self.reconstruction))
ax.set_title('Reconstruction Phase')
ax = fig.add_subplot(2,2,3 )
ax.imshow(np.abs(self.aperture))
ax.set_title("Aperture Magnitude")
ax = fig.add_subplot(2,2,4 )
ax.imshow(np.angle(self.aperture))
ax.set_title("Aperture Phase")
fig.tight_layout() Example 11
def phase_to_color_wheel(complex_number):
"""Map a phase of a complexnumber to a color in (r,g,b).
complex_number is phase is first mapped to angle in the range
[0, 2pi] and then to a color wheel with blue at zero phase.
"""
angles = np.angle(complex_number)
angle_round = int(((angles + 2 * np.pi) % (2 * np.pi))/np.pi*6)
# print(angleround)
color_map = {0: (0, 0, 1), # blue,
1: (0.5, 0, 1), # blue-violet
2: (1, 0, 1), # violet
3: (1, 0, 0.5), # red-violet,
4: (1, 0, 0), # red
5: (1, 0.5, 0), # red-oranage,
6: (1, 1, 0), # orange
7: (0.5, 1, 0), # orange-yellow
8: (0, 1, 0), # yellow,
9: (0, 1, 0.5), # yellow-green,
10: (0, 1, 1), # green,
11: (0, 0.5, 1) # green-blue,
}
return color_map[angle_round] Example 12
def time_sync(self, signal, start_est, stop_est):
self.symbol_found = 0
self.symbol_start = 0
self.log.debug('time_sync()')
corr_res = np.zeros( self.sym_len[self.mode], dtype=complex)
for s in range(start_est, self.sym_len[self.mode]-1):
corr_res[s] = self.gi_correlation(self.gi_len[self.mode], self.sym_len[self.mode], signal[s: s+self.sym_len[self.mode]] )
corr_max = np.argmax( np.abs(corr_res) )
self.symbol_found = 1
self.symbol_start = corr_max + 256 - 20
self.generated_samples = self.sym_len[self.mode] - self.gi_len[self.mode]
# Tracking
self.fine_freq_off = np.angle(corr_res[corr_max])/(2*np.pi*(self.sym_len[self.mode] - self.gi_len[self.mode])*(1.0/self.fs))
return Example 13
def compute_by_noise_pow(self, signal, n_pow):
s_spec = np.fft.fftpack.fft(signal * self._window)
s_amp = np.absolute(s_spec)
s_phase = np.angle(s_spec)
gamma = self._calc_aposteriori_snr(s_amp, n_pow)
xi = self._calc_apriori_snr(gamma)
self._prevGamma = gamma
nu = gamma * xi / (1.0 + xi)
self._G = (self._gamma15 * np.sqrt(nu) / gamma) * np.exp(-nu / 2.0) *\
((1.0 + nu) * spc.i0(nu / 2.0) + nu * spc.i1(nu / 2.0))
idx = np.less(s_amp ** 2.0, n_pow)
self._G[idx] = self._constant
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = xi[idx] / (xi[idx] + 1.0)
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = self._constant
self._G = np.maximum(self._G, 0.0)
amp = self._G * s_amp
amp = np.maximum(amp, 0.0)
amp2 = self._ratio * amp + (1.0 - self._ratio) * s_amp
self._prevAmp = amp
spec = amp2 * np.exp(s_phase * 1j)
return np.real(np.fft.fftpack.ifft(spec)) Example 14
def compute_by_noise_pow(self, signal, n_pow):
s_spec = np.fft.fftpack.fft(signal * self._window)
s_amp = np.absolute(s_spec)
s_phase = np.angle(s_spec)
gamma = self._calc_aposteriori_snr(s_amp, n_pow)
xi = self._calc_apriori_snr(gamma)
# xi = self._calc_apriori_snr2(gamma,n_pow)
self._prevGamma = gamma
nu = gamma * xi / (1.0 + xi)
self._G = xi / (1.0 + xi) * np.exp(0.5 * spc.exp1(nu))
idx = np.less(s_amp ** 2.0, n_pow)
self._G[idx] = self._constant
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = xi[idx] / (xi[idx] + 1.0)
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = self._constant
self._G = np.maximum(self._G, 0.0)
amp = self._G * s_amp
amp = np.maximum(amp, 0.0)
amp2 = self._ratio * amp + (1.0 - self._ratio) * s_amp
self._prevAmp = amp
spec = amp2 * np.exp(s_phase * 1j)
return np.real(np.fft.fftpack.ifft(spec)) Example 15
def compute_by_noise_pow(self, signal, n_pow):
s_spec = np.fft.fftpack.fft(signal * self._window)
s_amp = np.absolute(s_spec)
s_phase = np.angle(s_spec)
gamma = self._calc_aposteriori_snr(s_amp, n_pow)
# xi = self._calc_apriori_snr2(gamma,n_pow)
xi = self._calc_apriori_snr(gamma)
self._prevGamma = gamma
u = 0.5 - self._mu / (4.0 * np.sqrt(gamma * xi))
self._G = u + np.sqrt(u ** 2.0 + self._tau / (gamma * 2.0))
idx = np.less(s_amp ** 2.0, n_pow)
self._G[idx] = self._constant
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = xi[idx] / (xi[idx] + 1.0)
idx = np.isnan(self._G) + np.isinf(self._G)
self._G[idx] = self._constant
self._G = np.maximum(self._G, 0.0)
amp = self._G * s_amp
amp = np.maximum(amp, 0.0)
amp2 = self._ratio * amp + (1.0 - self._ratio) * s_amp
self._prevAmp = amp
spec = amp2 * np.exp(s_phase * 1j)
return np.real(np.fft.fftpack.ifft(spec)) Example 16
def plot_eigen_function(ax, se, n, psi1l, psi1r):
# plt.figure(figsize=(8, 8))
for x in range(se.info['L']):
for y in range(se.info['W']):
i = x + y * se.info['L']
w = np.sqrt(10) * np.abs(psi1r[i, n])
arg = np.angle(psi1r[i, n])
circle = plt.Circle((x, y), w, color='black', zorder=10)
ax.add_artist(circle)
ax.plot([x, x + w * np.cos(arg)], [y, y + w * np.sin(arg)], color='white', lw=.8, zorder=12)
ax.set_xlim([-.5, se.info['L'] - .5])
ax.set_ylim([-.5, se.info['W'] - .5])
# plt.show() Example 17
def single_spectrogram(inseq,fs,wlen,h,imag=False):
"""
imag: Return Imaginary Data of the STFT on True
"""
NFFT = int(2**(np.ceil(np.log2(wlen))))
K = np.sum(hamming(wlen, False))/wlen
raw_data = inseq.astype('float32')
raw_data = raw_data/np.amax(np.absolute(raw_data))
stft_data,_,_ = STFT(raw_data,wlen,h,NFFT,fs)
s = np.absolute(stft_data)/wlen/K;
if np.fmod(NFFT,2):
s[1:,:] *=2
else:
s[1:-2] *=2
real_data = np.transpose(20*np.log10(s + 10**-6)).astype(np.float32)
if imag:
imag_data = np.angle(stft_data).astype(np.float32)
return real_data,imag_data
return real_data Example 18
def complex2rgbalin(s): """ Displays complex image with intensity corresponding to the MODULUS and color (hsv) correponging to PHASE. From: pyVincent/ptycho.py """ ph=np.angle(s) t=np.pi/3 nx,ny=s.shape rgba=np.zeros((nx,ny,4)) rgba[:,:,0]=(ph<t)*(ph>-t) + (ph>t)*(ph<2*t)*(2*t-ph)/t + (ph>-2*t)*(ph<-t)*(ph+2*t)/t rgba[:,:,1]=(ph>t) + (ph<-2*t) *(-2*t-ph)/t+ (ph>0)*(ph<t) *ph/t rgba[:,:,2]=(ph<-t) + (ph>-t)*(ph<0) *(-ph)/t + (ph>2*t) *(ph-2*t)/t a=np.abs(s) a/=a.max() rgba[:,:,3]=a return rgba
Example 19
def autocorr(trace):
"""This function takes an obspy trace object and performs a phase autocorrelation
of the trace with itself. First, the hilbert transform is taken to obtain the
analytic signal and hence the instantaneous phase. This is then passed to a
fortran script where phase correlation is performed after Schimmel et al., 2011.
This function relies on the shared object file phasecorr.so, which is the file
containing the fortran subroutine for phase correlation.
"""
import numpy as np
from scipy.signal import hilbert
from phasecorr import phasecorr
# Hilbert transform to obtain the analytic signal
htrans = hilbert(trace)
# Perform phase autocorrelation with instantaneous phase
pac = phasecorr(np.angle(htrans),np.angle(htrans),len(htrans))
return pac Example 20
def __init__(self, gain):
self.gain = gain * np.pi # pi term puts angle output to pi range
# components
self.conjugate = Conjugate()
self.complex_mult = ComplexMultiply()
self.angle = Angle()
self.out = Sfix()
# specify component delay
self._delay = self.conjugate._delay + \
self.complex_mult._delay + \
self.angle._delay + 1
# constants
self.gain_sfix = Const(Sfix(self.gain, 3, -14)) Example 21
def NMF(stft, n_sources):
"""
Sound source separation using NMF
:param stft: the short-time Fourier Transform of the signal
:param n_sources: the number of sources
:returns:
- Ys: sources
"""
print "Computing approximations"
W, H = librosa.decompose.decompose(np.abs(stft), n_components=n_sources, sort=True)
print "Reconstructing signals"
Ys = list(scipy.outer(W[:,i], H[i])*np.exp(1j*np.angle(stft)) for i in xrange(n_sources))
print "Saving sound arrays"
ys = [librosa.istft(Y) for Y in Ys]
del Ys
return ys Example 22
def undo_melspect(spect, sample_rate=SAMPLE_RATE, fps=FPS, frame_len=FRAME_LEN, min_freq=MIN_FREQ, max_freq=MAX_FREQ, invert_melbank_method='transpose', phases='random', normalize=False):
"""
Resynthesizes a mel spectrogram into a numpy array of samples.
"""
# undo logarithmic scaling
spect = undo_logarithmize(spect)
# undo Mel filterbank
spect = undo_melfilter(spect, sample_rate, frame_len, min_freq, max_freq, invert_melbank_method)
# randomize or reuse phases
if phases == 'random':
spect = spect * np.exp(np.pi*2.j*np.random.random(spect.shape))
elif phases is not None:
spect = spect * np.exp(1.j * np.angle(phases))
# undo STFT
hop_size = sample_rate / fps
samples = undo_stft(spect, hop_size, frame_len)
# normalize if needed
if normalize:
samples -= samples.mean()
samples /= np.abs(samples).max()
return samples.astype(np.float32) Example 23
def __new__(cls, realpart, imagpart=None):
"""Create a new EMArray."""
# Create complex obj
if np.any(imagpart):
obj = np.real(realpart) + 1j*np.real(imagpart)
else:
obj = np.asarray(realpart, dtype=complex)
# Ensure its at least a 1D-Array, view cls
obj = np.atleast_1d(obj).view(cls)
# Store amplitude
obj.amp = np.abs(obj)
# Calculate phase, unwrap it, transform to degrees
obj.pha = np.rad2deg(np.unwrap(np.angle(obj.real + 1j*obj.imag)))
return obj
# 2. Input parameter checks for modelling
# 2.a <Check>s (alphabetically) Example 24
def ph_dec(m_phs, m_phc, mode='angle'):
if mode == 'sign':
m_bs = np.arcsin(m_phs)
m_bc = np.arccos(m_phc)
m_ph = np.sign(m_bs) * np.abs(m_bc)
elif mode == 'angle':
m_ph = np.angle(m_phc + m_phs * 1j)
return m_ph
#==============================================================================
# From 'analysis_with_del_comp_and_ph_encoding_from_files'
# f0_type: 'f0', 'lf0' Example 25
def phase_string(sig):
"""Take the angle and create a string as \pi if close to some \pi values"""
if isinstance(sig, np.ndarray):
sig = sig.ravel()[0]
if isinstance(sig, np.complex):
angle = np.angle(sig)
else:
angle = sig
fractions = [Fraction(i, 12) for i in np.arange(-12, 12 + 1)]
pi_multiples = np.array([np.pi * frac_to_float(f) for f in fractions])
strings = [frac_to_str(f) for f in fractions]
if np.min(np.abs(angle - pi_multiples)) < 1e-3:
return strings[np.argmin(np.abs(angle - pi_multiples))]
else:
return '%.2f' % angle Example 26
def cimshow(f_impulse):
plt.figure(figsize=(7,5))
plt.subplot(2,2,1)
plt.imshow(np.real(f_impulse))
plt.colorbar()
plt.title('Re{}')
plt.subplot(2,2,2)
plt.imshow(np.imag(f_impulse))
plt.colorbar()
plt.title('Img{}')
plt.subplot(2,2,2+1)
plt.imshow(np.abs(f_impulse))
plt.colorbar()
plt.title('Magnitude')
plt.subplot(2,2,2+2)
plt.imshow(np.angle(f_impulse))
plt.colorbar()
plt.title('Phase') Example 27
def cimshow(f_impulse):
plt.figure(figsize=(7,5))
plt.subplot(2,2,1)
plt.imshow(np.real(f_impulse))
plt.colorbar()
plt.title('Re{}')
plt.subplot(2,2,2)
plt.imshow(np.imag(f_impulse))
plt.colorbar()
plt.title('Img{}')
plt.subplot(2,2,2+1)
plt.imshow(np.abs(f_impulse))
plt.colorbar()
plt.title('Magnitude')
plt.subplot(2,2,2+2)
plt.imshow(np.angle(f_impulse))
plt.colorbar()
plt.title('Phase') Example 28
def cimshow(f_impulse):
plt.figure(figsize=(7,5))
plt.subplot(2,2,1)
plt.imshow(np.real(f_impulse))
plt.colorbar()
plt.title('Re{}')
plt.subplot(2,2,2)
plt.imshow(np.imag(f_impulse))
plt.colorbar()
plt.title('Img{}')
plt.subplot(2,2,2+1)
plt.imshow(np.abs(f_impulse))
plt.colorbar()
plt.title('Magnitude')
plt.subplot(2,2,2+2)
plt.imshow(np.angle(f_impulse))
plt.colorbar()
plt.title('Phase') Example 29
def update_recon_py(Recon1, support):
err1 = 1.0
Constraint = np.ones(Recon1.shape)
# cdef int R1y, R1x
Recon1_abs = np.abs(Recon1)
Recon1_pwr2 = np.power(Recon1_abs, 2)
R1y, R1x = Recon1.shape
for p in range(R1y):
for q in range(R1x):
if support[p, q] == 1:
Constraint[p, q] = Recon1_abs[p, q]
err1 += Recon1_pwr2[p, q]
if Recon1_abs[p, q] > 1:
Constraint[p, q] = 1
Recon1_update = Constraint * np.exp(1j * np.angle(Recon1))
return Recon1_update, err1 Example 30
def imshow_GfpGbp(Gfp, Gbp):
plt.subplot(2, 2, 1)
plt.imshow(np.abs(Gfp), cmap='Greys_r')
plt.title('abs(Gfp)')
plt.subplot(2, 2, 2)
plt.imshow(np.angle(Gfp), cmap='Greys_r')
plt.title('angle(Gfp)')
plt.subplot(2, 2, 1 + 2)
plt.imshow(np.abs(Gbp), cmap='Greys_r')
plt.title('abs(Gbp)')
plt.subplot(2, 2, 2 + 2)
plt.imshow(np.angle(Gbp), cmap='Greys_r')
plt.title('angle(Gbp)') Example 31
def imshow_h(self):
Gbp, Gfp = self.Gbp, self.Gfp
plt.figure(figsize=[10, 10])
plt.subplot(2, 2, 1)
plt.imshow(np.abs(Gbp), cmap='Greys_r', interpolation='none')
plt.title('Gbp - Maganitute')
plt.subplot(2, 2, 2)
plt.imshow(np.angle(Gbp), cmap='Greys_r', interpolation='none')
plt.title('Gbp - Angle')
plt.subplot(2, 2, 2 + 1)
plt.imshow(np.abs(Gfp), cmap='Greys_r', interpolation='none')
plt.title('Gbf - Maganitute')
plt.subplot(2, 2, 2 + 2)
plt.imshow(np.angle(Gfp), cmap='Greys_r', interpolation='none')
plt.title('Gbf - Angle') Example 32
def run_complex(self, Original):
"""
diffract original image and recon.
ALso, the results will be shown
"""
Input = self.diffract_full(Original)
Recon = self.reconstruct(Input)
figure(figsize=(3 * 3 + 2, 3))
subplot(1, 4, 1)
imshow(Original, 'Greys_r')
title('Original')
subplot(1, 4, 2)
imshow(np.abs(Input), 'Greys_r')
title('|Diffraction|')
subplot(1, 4, 3)
imshow(np.angle(Input), 'Greys_r')
title('Angle(Diffraction)')
subplot(1, 4, 4)
imshow(np.abs(Recon), 'Greys_r')
title('Modulus: |Recon|') Example 33
def fourierExtrapolation(x, n_predict):
n = len(x)
n_harm = 10 # number of harmonics in model
t = np.arange(0, n)
p = np.polyfit(t, x, 1) # find linear trend in x
x_notrend = x - p[0] * t # detrended x
x_freqdom = fft.fft(x_notrend) # detrended x in frequency domain
f = fft.fftfreq(n) # frequencies
indexes = list(range(n))
# sort indexes by frequency, lower -> higher
indexes.sort(key = lambda i: np.absolute(f[i]))
t = np.arange(0, n + n_predict)
restored_sig = np.zeros(t.size)
for i in indexes[:1 + n_harm * 2]:
ampli = np.absolute(x_freqdom[i]) / n # amplitude
phase = np.angle(x_freqdom[i]) # phase
restored_sig += ampli * np.cos(2 * np.pi * f[i] * t + phase)
return restored_sig + p[0] * t Example 34
def fourierExtrapolation(x, n_predict):
n = len(x)
n_harm = 10 # number of harmonics in model
t = np.arange(0, n)
p = np.polyfit(t, x, 1) # find linear trend in x
x_notrend = x - p[0] * t # detrended x
x_freqdom = fft.fft(x_notrend) # detrended x in frequency domain
f = fft.fftfreq(n) # frequencies
indexes = list(range(n))
# sort indexes by frequency, lower -> higher
indexes.sort(key = lambda i: np.absolute(f[i]))
t = np.arange(0, n + n_predict)
restored_sig = np.zeros(t.size)
for i in indexes[:1 + n_harm * 2]:
ampli = np.absolute(x_freqdom[i]) / n # amplitude
phase = np.angle(x_freqdom[i]) # phase
restored_sig += ampli * np.cos(2 * np.pi * f[i] * t + phase)
return restored_sig + p[0] * t Example 35
def fourierExtrapolation(x, n_predict):
n = len(x)
n_harm = 10 # number of harmonics in model
t = np.arange(0, n)
p = np.polyfit(t, x, 1) # find linear trend in x
x_notrend = x - p[0] * t # detrended x
x_freqdom = fft.fft(x_notrend) # detrended x in frequency domain
f = fft.fftfreq(n) # frequencies
indexes = list(range(n))
# sort indexes by frequency, lower -> higher
indexes.sort(key = lambda i: np.absolute(f[i]))
t = np.arange(0, n + n_predict)
restored_sig = np.zeros(t.size)
for i in indexes[:1 + n_harm * 2]:
ampli = np.absolute(x_freqdom[i]) / n # amplitude
phase = np.angle(x_freqdom[i]) # phase
restored_sig += ampli * np.cos(2 * np.pi * f[i] * t + phase)
return restored_sig + p[0] * t Example 36
def griffinlim(spectrogram, n_iter=50, window='hann', n_fft=2048, win_length=2048, hop_length=-1, verbose=False):
if hop_length == -1:
hop_length = n_fft // 4
angles = np.exp(2j * np.pi * np.random.rand(*spectrogram.shape))
t = tqdm(range(n_iter), ncols=100, mininterval=2.0, disable=not verbose)
for i in t:
full = np.abs(spectrogram).astype(np.complex) * angles
inverse = librosa.istft(full, hop_length = hop_length, win_length = win_length, window = window)
rebuilt = librosa.stft(inverse, n_fft = n_fft, hop_length = hop_length, win_length = win_length, window = window)
angles = np.exp(1j * np.angle(rebuilt))
if verbose:
diff = np.abs(spectrogram) - np.abs(rebuilt)
t.set_postfix(loss=np.linalg.norm(diff, 'fro'))
full = np.abs(spectrogram).astype(np.complex) * angles
inverse = librosa.istft(full, hop_length = hop_length, win_length = win_length, window = window)
return inverse Example 37
def fourierEx(x, n_predict, harmonics):
n = len(x) # number of input samples
n_harmonics = harmonics # f, 2*f, 3*f, .... n_harmonics ( 1,2, )
t = np.arange(0, n) # place to store data
p = np.polyfit(t, x, 1) # find trend
x_no_trend = x - p[0] * t
x_frequency_domains = fft.fft(x_no_trend)
f = np.fft.fftfreq(n) # frequencies
indexes = list(range(n))
indexes.sort(key=lambda i: np.absolute(f[i]))
t = np.arange(0, n + n_predict)
restored_signal = np.zeros(t.size)
for i in indexes[:1 + n_harmonics * 2]:
amplitude = np.absolute(x_frequency_domains[i] / n)
phase = np.angle(x_frequency_domains[i])
restored_signal += amplitude * np.cos(2 * np.pi * f[i] * t + phase)
return restored_signal + p[0] * t Example 38
def processData(self, data):
times = data.xvals('Time')
dt = times[1]-times[0]
data1 = data.asarray()
ft = np.fft.fft(data1)
## determine frequencies in fft data
df = 1.0 / (len(data1) * dt)
freqs = np.linspace(0.0, (len(ft)-1) * df, len(ft))
## flatten spikes at f0 and harmonics
f0 = self.ctrls['f0'].value()
for i in xrange(1, self.ctrls['harmonics'].value()+2):
f = f0 * i # target frequency
## determine index range to check for this frequency
ind1 = int(np.floor(f / df))
ind2 = int(np.ceil(f / df)) + (self.ctrls['samples'].value()-1)
if ind1 > len(ft)/2.:
break
mag = (abs(ft[ind1-1]) + abs(ft[ind2+1])) * 0.5
for j in range(ind1, ind2+1):
phase = np.angle(ft[j]) ## Must preserve the phase of each point, otherwise any transients in the trace might lead to large artifacts.
re = mag * np.cos(phase)
im = mag * np.sin(phase)
ft[j] = re + im*1j
ft[len(ft)-j] = re - im*1j
data2 = np.fft.ifft(ft).real
ma = metaarray.MetaArray(data2, info=data.infoCopy())
return ma Example 39
def removePeriodic(data, f0=60.0, dt=None, harmonics=10, samples=4):
if (hasattr(data, 'implements') and data.implements('MetaArray')):
data1 = data.asarray()
if dt is None:
times = data.xvals('Time')
dt = times[1]-times[0]
else:
data1 = data
if dt is None:
raise Exception('Must specify dt for this data')
ft = np.fft.fft(data1)
## determine frequencies in fft data
df = 1.0 / (len(data1) * dt)
freqs = np.linspace(0.0, (len(ft)-1) * df, len(ft))
## flatten spikes at f0 and harmonics
for i in xrange(1, harmonics + 2):
f = f0 * i # target frequency
## determine index range to check for this frequency
ind1 = int(np.floor(f / df))
ind2 = int(np.ceil(f / df)) + (samples-1)
if ind1 > len(ft)/2.:
break
mag = (abs(ft[ind1-1]) + abs(ft[ind2+1])) * 0.5
for j in range(ind1, ind2+1):
phase = np.angle(ft[j]) ## Must preserve the phase of each point, otherwise any transients in the trace might lead to large artifacts.
re = mag * np.cos(phase)
im = mag * np.sin(phase)
ft[j] = re + im*1j
ft[len(ft)-j] = re - im*1j
data2 = np.fft.ifft(ft).real
if (hasattr(data, 'implements') and data.implements('MetaArray')):
return metaarray.MetaArray(data2, info=data.infoCopy())
else:
return data2 Example 40
def processData(self, data):
times = data.xvals('Time')
dt = times[1]-times[0]
data1 = data.asarray()
ft = np.fft.fft(data1)
## determine frequencies in fft data
df = 1.0 / (len(data1) * dt)
freqs = np.linspace(0.0, (len(ft)-1) * df, len(ft))
## flatten spikes at f0 and harmonics
f0 = self.ctrls['f0'].value()
for i in xrange(1, self.ctrls['harmonics'].value()+2):
f = f0 * i # target frequency
## determine index range to check for this frequency
ind1 = int(np.floor(f / df))
ind2 = int(np.ceil(f / df)) + (self.ctrls['samples'].value()-1)
if ind1 > len(ft)/2.:
break
mag = (abs(ft[ind1-1]) + abs(ft[ind2+1])) * 0.5
for j in range(ind1, ind2+1):
phase = np.angle(ft[j]) ## Must preserve the phase of each point, otherwise any transients in the trace might lead to large artifacts.
re = mag * np.cos(phase)
im = mag * np.sin(phase)
ft[j] = re + im*1j
ft[len(ft)-j] = re - im*1j
data2 = np.fft.ifft(ft).real
ma = metaarray.MetaArray(data2, info=data.infoCopy())
return ma Example 41
def removePeriodic(data, f0=60.0, dt=None, harmonics=10, samples=4):
if (hasattr(data, 'implements') and data.implements('MetaArray')):
data1 = data.asarray()
if dt is None:
times = data.xvals('Time')
dt = times[1]-times[0]
else:
data1 = data
if dt is None:
raise Exception('Must specify dt for this data')
ft = np.fft.fft(data1)
## determine frequencies in fft data
df = 1.0 / (len(data1) * dt)
freqs = np.linspace(0.0, (len(ft)-1) * df, len(ft))
## flatten spikes at f0 and harmonics
for i in xrange(1, harmonics + 2):
f = f0 * i # target frequency
## determine index range to check for this frequency
ind1 = int(np.floor(f / df))
ind2 = int(np.ceil(f / df)) + (samples-1)
if ind1 > len(ft)/2.:
break
mag = (abs(ft[ind1-1]) + abs(ft[ind2+1])) * 0.5
for j in range(ind1, ind2+1):
phase = np.angle(ft[j]) ## Must preserve the phase of each point, otherwise any transients in the trace might lead to large artifacts.
re = mag * np.cos(phase)
im = mag * np.sin(phase)
ft[j] = re + im*1j
ft[len(ft)-j] = re - im*1j
data2 = np.fft.ifft(ft).real
if (hasattr(data, 'implements') and data.implements('MetaArray')):
return metaarray.MetaArray(data2, info=data.infoCopy())
else:
return data2 Example 42
def getMag(_fft):
_mag = abs(_fft);
norm_mag = _mag/_mag.max();
#angle = angle/angle.max();
return norm_mag; Example 43
def resample_2d_complex(array, sample_pts, query_pts):
''' Resamples a 2D complex array by interpolating the magnitude and phase
independently and merging the results into a complex value.
Args:
array (numpy.ndarray): complex 2D array.
sample_pts (tuple): pair of numpy.ndarray objects that contain the x and y sample locations,
each array should be 1D.
query_pts (tuple): points to interpolate onto, also 1D for each array.
Returns:
numpy.ndarray array resampled onto query_pts via bivariate spline
'''
xq, yq = np.meshgrid(*query_pts)
mag = abs(array)
phase = np.angle(array)
magfunc = interpolate.RegularGridInterpolator(sample_pts, mag)
phasefunc = interpolate.RegularGridInterpolator(sample_pts, phase)
interp_mag = magfunc((yq, xq))
interp_phase = phasefunc((yq, xq))
return interp_mag * exp(1j * interp_phase) Example 44
def phase(self):
"""
Return the phase spectrum.
"""
return np.angle(self) Example 45
def phase(z):
val = np.angle(z)
# val = np.rad2deg(np.unwrap(np.angle((z))))
return val Example 46
def DFT(x, w, N):
""" Discrete Fourier Transformation(Analysis) of a given real input signal
via an FFT implementation from scipy. Single channel is being supported.
Args:
x : (array) Real time domain input signal
w : (array) Desired windowing function
N : (int) FFT size
Returns:
magX : (2D ndarray) Magnitude Spectrum
phsX : (2D ndarray) Phase Spectrum
"""
# Half spectrum size containing DC component
hlfN = (N/2)+1
# Half window size. Two parameters to perform zero-phase windowing technique
hw1 = int(math.floor((w.size+1)/2))
hw2 = int(math.floor(w.size/2))
# Window the input signal
winx = x*w
# Initialize FFT buffer with zeros and perform zero-phase windowing
fftbuffer = np.zeros(N)
fftbuffer[:hw1] = winx[hw2:]
fftbuffer[-hw2:] = winx[:hw2]
# Compute DFT via scipy's FFT implementation
X = fft(fftbuffer)
# Acquire magnitude and phase spectrum
magX = (np.abs(X[:hlfN]))
phsX = (np.angle(X[:hlfN]))
return magX, phsX Example 47
def test_simple(self):
x = np.array([1+0j, 1+2j])
y_r = np.log(np.abs(x)) + 1j * np.angle(x)
y = np.log(x)
for i in range(len(x)):
assert_almost_equal(y[i], y_r[i]) Example 48
def test_basic_ufuncs(self):
# Test various functions such as sin, cos.
(x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
assert_equal(np.cos(x), cos(xm))
assert_equal(np.cosh(x), cosh(xm))
assert_equal(np.sin(x), sin(xm))
assert_equal(np.sinh(x), sinh(xm))
assert_equal(np.tan(x), tan(xm))
assert_equal(np.tanh(x), tanh(xm))
assert_equal(np.sqrt(abs(x)), sqrt(xm))
assert_equal(np.log(abs(x)), log(xm))
assert_equal(np.log10(abs(x)), log10(xm))
assert_equal(np.exp(x), exp(xm))
assert_equal(np.arcsin(z), arcsin(zm))
assert_equal(np.arccos(z), arccos(zm))
assert_equal(np.arctan(z), arctan(zm))
assert_equal(np.arctan2(x, y), arctan2(xm, ym))
assert_equal(np.absolute(x), absolute(xm))
assert_equal(np.angle(x + 1j*y), angle(xm + 1j*ym))
assert_equal(np.angle(x + 1j*y, deg=True), angle(xm + 1j*ym, deg=True))
assert_equal(np.equal(x, y), equal(xm, ym))
assert_equal(np.not_equal(x, y), not_equal(xm, ym))
assert_equal(np.less(x, y), less(xm, ym))
assert_equal(np.greater(x, y), greater(xm, ym))
assert_equal(np.less_equal(x, y), less_equal(xm, ym))
assert_equal(np.greater_equal(x, y), greater_equal(xm, ym))
assert_equal(np.conjugate(x), conjugate(xm)) Example 49
def synthesis_speech(noisy_speech, ideal_mask, win_type, win_len, shift_len, syn_method='A&R'):
samples = noisy_speech.shape[0]
frames = (samples - win_len) // shift_len
if win_type == 'hanning':
window = np.hanning(win_len)
elif win_type == 'hamming':
window = np.hamming(win_len)
elif win_type == 'rectangle':
window = np.ones(win_len)
to_ifft = np.zeros(win_len, dtype=np.complex64)
clean_speech = np.zeros((frames-1)*shift_len+win_len, dtype=np.float32)
window_sum = np.zeros((frames-1)*shift_len+win_len, dtype=np.float32)
for i in range(frames):
one_frame = noisy_speech[i * shift_len: i * shift_len + win_len]
windowed_frame = np.multiply(one_frame, window)
stft = np.fft.fft(windowed_frame, win_len)
masked_abs = np.abs(stft[:win_len//2+1]) * ideal_mask[:, i]
to_ifft[:win_len//2+1] = masked_abs * np.exp(1j * np.angle(stft[:win_len//2+1]))
to_ifft[win_len//2+1:] = np.conj(to_ifft[win_len//2-1:0:-1])
speech_seg = np.real(np.fft.ifft(to_ifft, win_len))
if syn_method == 'A&R' or syn_method == 'ALLEN & RABINER':
clean_speech[i*shift_len:i*shift_len+win_len] += speech_seg
window_sum[i*shift_len:i*shift_len+win_len] += window
elif syn_method == 'G&L' or syn_method == 'GRIFFIN & LIM':
speech_seg = np.multiply(speech_seg, window)
clean_speech[i * shift_len:i * shift_len + win_len] += speech_seg
window_sum[i * shift_len:i * shift_len + win_len] += np.power(window, 2.)
# if i > 0:
# clean_speech[i*shift_len: (i-1)*shift_len+win_len] *= 0.5
window_sum = np.where(window_sum < 1e-2, 1e-2, window_sum)
return clean_speech / window_sum Example 50
def synthesis_speech(ns, mk, win_type, win_len, shift_len, syn_method='A&R'):
samples = ns.shape[0]
frames = (samples - win_len) // shift_len
if win_type == 'hanning':
window = np.hanning(win_len)
elif win_type == 'hamming':
window = np.hamming(win_len)
elif win_type == 'rectangle':
window = np.ones(win_len)
to_ifft = np.zeros(win_len, dtype=np.complex64)
clean_speech = np.zeros((frames-1)*shift_len+win_len, dtype=np.float32)
window_sum = np.zeros((frames-1)*shift_len+win_len, dtype=np.float32)
for i in range(frames):
one_frame = ns[i * shift_len: i * shift_len + win_len]
windowed_frame = np.multiply(one_frame, window)
stft = np.fft.fft(windowed_frame, win_len)
masked_abs = np.abs(stft[:win_len//2+1]) * mk[:, i]
to_ifft[:win_len//2+1] = masked_abs * np.exp(1j * np.angle(stft[:win_len//2+1]))
to_ifft[win_len//2+1:] = np.conj(to_ifft[win_len//2-1:0:-1])
speech_seg = np.real(np.fft.ifft(to_ifft, 320))
if syn_method == 'A&R' or syn_method == 'ALLEN & RABINER':
clean_speech[i*shift_len:i*shift_len+win_len] += speech_seg
window_sum[i*shift_len:i*shift_len+win_len] += window
elif syn_method == 'G&L' or syn_method == 'GRIFFIN & LIM':
speech_seg = np.multiply(speech_seg, window)
clean_speech[i * shift_len:i * shift_len + win_len] += speech_seg
window_sum[i * shift_len:i * shift_len + win_len] += np.power(window, 2.)
# if i > 0:
# clean_speech[i*shift_len: (i-1)*shift_len+win_len] *= 0.5
window_sum = np.where(window_sum < 1e-2, 1e-2, window_sum)
return clean_speech / window_sum 