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 _inside_contour(pos, contour):
"""Aux function"""
npos = len(pos)
x, y = pos[:, :2].T
check_mask = np.ones((npos), dtype=bool)
check_mask[((x < np.min(x)) | (y < np.min(y)) |
(x > np.max(x)) | (y > np.max(y)))] = False
critval = 0.1
sel = np.where(check_mask)[0]
for this_sel in sel:
contourx = contour[:, 0] - pos[this_sel, 0]
contoury = contour[:, 1] - pos[this_sel, 1]
angle = np.arctan2(contoury, contourx)
angle = np.unwrap(angle)
total = np.sum(np.diff(angle))
check_mask[this_sel] = np.abs(total) > critval
return check_mask Example 2
def plot_poses(poses, pose_type='absolute'):
f = plt.figure(1)
f.clf()
f.subplots_adjust(hspace=0.25)
# , **{'ylabel': 'Position X(m)', 'xlim': [min(xs), max(xs)]}
pparams = { 'linewidth':1, }
# font = { 'family': 'normal', 'weight': 'normal', 'size': 12 }
# mpl.rc('font', **font)
rpyxyz = np.vstack([pose.to_rpyxyz() for pose in poses])
ax = f.add_subplot(2, 1, 1)
if pose_type == 'absolute':
for j,label in enumerate(['Roll', 'Pitch', 'Yaw']):
ax.plot(np.unwrap(rpyxyz[:,j], discont=np.pi), label=label)
elif pose_type == 'relative':
for j,label in enumerate(['Roll', 'Pitch', 'Yaw']):
ax.plot(np.unwrap(rpyxyz[1:,j]-rpyxyz[:-1,j], discont=np.pi), label=label)
else:
raise RuntimeError('Unknown pose_type=%s, use absolute or relative' % pose_type)
ax.legend(loc='upper right', fancybox=False, ncol=3,
# bbox_to_anchor=(1.0,1.2),
prop={'size':13})
ax = f.add_subplot(2, 1, 2)
for j,label in enumerate(['X', 'Y', 'Z']):
ax.plot(rpyxyz[:,j+3], label=label)
ax.legend(loc='upper right', fancybox=False, ncol=3,
# bbox_to_anchor=(1.0,1.2),
prop={'size':13})
plt.tight_layout()
plt.show(block=True) Example 3
def another():
# plot the figure of signals
fig = plt.figure(1)
ax = fig.add_subplot(311)
ax.set_title('input signal')
ax.plot(t[1:Npts],x[1:Npts]) # just plot part of the signal
ax = fig.add_subplot(312)
ax.set_title('expected signal')
ax.plot(t[1:Npts],xo[1:Npts])
ax = fig.add_subplot(313)
ax.set_title('filtered signal')
ax.plot(t[1:Npts],y[1:Npts])
fig.savefig('pic1.png')
# plot the mag & phase response of the LPF
w,h = signal.freqz(b,1)
hdb = 20 * np.log10(np.abs(h))
hphs = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
fig2=plt.figure(2)
ax2 = fig2.add_subplot(211)
ax2.set_title('frequency response')
ax2.plot(w/max(w),hdb)
ax2 = fig2.add_subplot(212)
ax2.set_title('phase response')
ax2.plot(w/max(w),hphs)
fig2.savefig('pic2.png') Example 4
def unwrap(plotar, ms2mappings):
for k in plotar.keys():
for d in plotar[k].keys():
# get a reference to the data set
dsref = plotar[k][d]
# get sorted data and unwrap the phases
(xvals, yvals) = zip(*sorted(zip(dsref.xval, dsref.yval), key=operator.itemgetter(0)))
yvals = numpy.unwrap(numpy.deg2rad(yvals))
dsref.xval = numpy.array(xvals)
dsref.yval = numpy.rad2deg(yvals) Example 5
def phaserate(plotar, ms2mappings):
if plotar.plotType!='phatime':
raise RuntimeError, "phaserate() cannot run on plot type {0}".format( plotar.plotType )
spm = ms2mappings.spectralMap
# iterate over all plots and all data sets within the plots
for k in plotar.keys():
for d in plotar[k].keys():
# get a reference to the data set
dsref = plotar[k][d]
# get the full data set label - we have access to all the data set's properties (FQ, SB, POL etc)
n = plots.join_label(k, d)
# fit a line through the unwrapped phase
unw = numpy.unwrap(numpy.deg2rad(dsref.yval))
coeffs = numpy.polyfit(dsref.xval, unw, 1)
# evaluate the fitted polynomial at the x-loci
extray = numpy.polyval(coeffs, dsref.xval)
# here we could compute the reliability of the fit
diff = unw - extray
ss_tot = numpy.sum(numpy.square(unw - unw.mean()))
ss_res = numpy.sum(numpy.square(diff))
r_sq = 1.0 - ss_res/ss_tot
# compare std deviation and variance in the residuals after fit
std_r = numpy.std(diff)
var_r = numpy.var(diff)
f = spm.frequencyOfFREQ_SB(n.FQ, n.SB)
rate = coeffs[0]
if var_r<std_r:
print "{0}: {1:.8f} ps/s @ {2:5.4f}MHz [R2={3:.3f}]".format(n, rate/(2.0*numpy.pi*f*1.0e-12), f/1.0e6, r_sq )
# before plotting wrap back to -pi,pi and transform to degrees
dsref.extra = [ drawline_fn(n, dsref.xval, numpy.rad2deg(do_wrap(extray))) ] Example 6
def phasedbg(plotar, ms2mappings):
global cache
if plotar.plotType!='phatime':
raise RuntimeError, "phasedbg() cannot run on plot type {0}".format( plotar.plotType )
store = len(cache)==0
# iterate over all plots and all data sets within the plots
for k in plotar.keys():
for d in plotar[k].keys():
# get a reference to the data set
dsref = plotar[k][d]
# get the full data set label - we have access to all the data set's properties (FQ, SB, POL etc)
n = plots.join_label(k, d)
# fit a line through the unwrapped phase
unw = numpy.unwrap(numpy.deg2rad(dsref.yval))
#coeffs = numpy.polyfit(dsref.xval, unw, 1)
coeffs = numpy.polyfit(xrange(len(dsref.yval)), unw, 1)
# evaluate the fitted polynomial at the x-loci
extray = numpy.polyval(coeffs, dsref.xval)
# here we could compute the reliability of the fit
diff = unw - extray
# compare std deviation and variance in the residuals after fit
std_r = numpy.std(diff)
var_r = numpy.var(diff)
coeffs = numpy.rad2deg(coeffs)
if var_r<std_r:
# decide what to do
if store:
cache[n] = coeffs
else:
# check if current key exists in cache; if so
# do differencing
otherCoeffs = cache.get(n, None)
if otherCoeffs is None:
print "{0}: not found in cache".format( n )
else:
delta = otherCoeffs - coeffs
print "{0.BL} {0.SB} {0.P}: dRate={1:5.4f} dOff={2:4.1f}".format(n, delta[0], delta[1]+360.0 if delta[1]<0.0 else delta[1])
# before plotting wrap back to -pi,pi and transform to degrees
#dsref.extra = [ drawline_fn(n, dsref.xval, numpy.rad2deg(do_wrap(extray))) ]
if not store:
cache = {} Example 7
def freqz(b, a=1, fs=2.0, worN=None, whole=False):
"""Plot frequency response of a filter.
This is a convenience function to plot frequency response, and internally uses
:func:`scipy.signal.freqz` to estimate the response. For further details, see the
documentation for :func:`scipy.signal.freqz`.
:param b: numerator of a linear filter
:param a: denominator of a linear filter
:param fs: sampling rate in Hz (optional, normalized frequency if not specified)
:param worN: see :func:`scipy.signal.freqz`
:param whole: see :func:`scipy.signal.freqz`
:returns: (frequency vector, frequency response vector)
>>> import arlpy
>>> arlpy.signal.freqz([1,1,1,1,1], fs=120000);
>>> w, h = arlpy.signal.freqz([1,1,1,1,1], fs=120000)
"""
import matplotlib.pyplot as plt
w, h = _sig.freqz(b, a, worN, whole)
f = w*fs/(2*_np.pi)
fig = plt.figure()
ax1 = fig.add_subplot(111)
plt.plot(f, 20*_np.log10(abs(h)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [Hz]')
plt.grid()
ax1.twinx()
angles = _np.unwrap(_np.angle(h))
plt.plot(f, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.axis('tight')
plt.show()
return w, h Example 8
def plotPhaseResponse(self, freqs=None, scale='radians',
showCorners=True,
label='Frequency Response',
ax=None, **kwargs):
"""Plot the frequency response of the filter.
"""
if ax is None:
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
freqs, responses = self.frequencyResponse(freqs=freqs)
freqs = freqs * self.sampRate * 0.5 / np.pi
responseAngles = np.unwrap(np.angle(responses))
scale = scale.lower()
if scale == 'radians':
ax.set_ylabel('Phase (Radians)')
elif scale == 'cycles':
ax.set_ylabel('Phase (Cycles)')
responseAngles /= 2.0*np.pi
elif scale == 'degrees':
ax.set_ylabel('Phase (Degrees)')
responseAngles = responseAngles*180.0 / np.pi
else:
raise Exception('Invalid scale: ' + str(scale) + '.')
lines = ax.plot(freqs, responseAngles,
label=label, **kwargs)
result = {'ax': ax, 'lines': lines}
if showCorners:
cornerLines = ax.vlines((self.lowFreq,self.highFreq),
np.min(responseAngles), np.max(responseAngles),
color='violet', linestyle='--', label='Corners')
result['cornerLines'] = cornerLines
ax.set_xlabel('Frequency (Hz)')
return result Example 9
def _calc(self, data, ret_obj, **kwargs):
self._inst_als = _AlsCvxopt(**kwargs)
try:
shp = data.shape[0:-2]
total_num = _np.array(shp).prod()
counter = 1
for idx in _np.ndindex(shp):
print('Detrended iteration {} / {}'.format(counter, total_num))
ph = _np.unwrap(_np.angle(data[idx]))
if self.rng is None:
err_phase = self._inst_als.calculate(ph)
else:
err_phase = self._inst_als.calculate(ph[..., self.rng])
h = _np.zeros(err_phase.shape)
h += _hilbert(err_phase)
correction_factor = 1/_np.exp(h) * _np.exp(-1j*err_phase)
if self.rng is None:
ret_obj[idx] *= correction_factor
else:
ret_obj[idx][..., self.rng] *= correction_factor
counter += 1
except:
return False
else:
# print(self._inst_als.__dict__)
return True Example 10
def phase_resp(self, frequencies=None, unwrap=False):
"""Calculate the real phase response"""
w, h = self.complex_freq_resp(frequencies)
phase = np.angle(h, deg=False)
phase = np.unwrap(phase) if unwrap else phase
phase = np.rad2deg(phase)
freqs = rad2hz(w, self.fs)
return freqs, phase Example 11
def plot_mag_phase(self, filename=None, plotpoints=10000, unwrap=False):
"""Produce a plot with magnitude and phase response in the same
figure. The y-axis on the left side belongs to the magnitude
response. The y-axis on the right side belongs to the phase
response
"""
unused_freq, mag = self.magnitude_resp(plotpoints)
freq, pha = self.phase_resp(plotpoints, unwrap=unwrap)
fig = plt.figure(1)
ax_mag = fig.add_subplot(111)
ax_pha = ax_mag.twinx()
ax_mag.semilogx(freq, mag, label='magnitude', color='red', linestyle='-')
ax_pha.semilogx(freq, pha, label='phase', color='blue', linestyle='--')
ax_mag.grid(True)
ax_mag.set_xlim(10, self.fs/2)
#ax_mag.set_ylim(bottom=-80) # FIXME: ad proper padding
#ax_mag.margins(0.1)
ax_mag.set_title('Frequency response')
ax_mag.set_xlabel('Frequency [Hz]')
ax_mag.set_ylabel('Magnitude [dB]')
ax_pha.set_ylabel('Phase [deg]')
handles1, labels1 = ax_mag.get_legend_handles_labels()
handles2, labels2 = ax_pha.get_legend_handles_labels()
plt.legend(handles1 + handles2, labels1 + labels2, loc='best')
if filename is None:
plt.show()
else:
try:
plt.savefig(filename)
finally:
plt.close(1) Example 12
def phaseunwrap(array):
"""
unwraps the phase from a given complex array returning a signal with 0 average phase slope
"""
data = np.asarray(array)
phase = np.unwrap(np.angle(data))
avg = np.average(np.diff(phase))
data = [x*np.exp(-1j*avg*i) for i,x in enumerate(data)]
# print(np.average(np.diff(np.angle(data))))
return np.asarray(data) Example 13
def find_resonance(frec,S11,margin=31,doplots=False):
"""
Returns the resonance frequency from maximum of the derivative of the phase
"""
#Smooth data for initial guesses
sReS11 = np.array(smooth(S11.real,margin,3))
sImS11 = np.array(smooth(S11.imag,margin,3))
sS11 = np.array( [x+1j*y for x,y in zip(sReS11,sImS11) ] )
#Make smoothed phase vector removing 2pi jumps
sArgS11 = np.angle(sS11)
sArgS11 = np.unwrap(sArgS11)
sdiffang = np.diff(sArgS11)
#Get resonance index from maximum of the derivative of the phase
avgph = np.average(sdiffang)
errvec = [np.power(x-avgph,2.) for x in sdiffang]
ires = np.argmax(errvec[margin:-margin])+margin
f0i=frec[ires]
print("Max index: ",ires," Max frequency: ",f0i)
if doplots:
plt.title('Original signal (Re,Im)')
plt.plot(frec,S11.real)
plt.plot(frec,S11.imag)
plt.show()
plt.plot(np.angle(sS11))
plt.title('Smoothed Phase')
plt.axis('auto')
plt.show()
plt.plot(sdiffang[margin:-margin])
plt.plot(sdiffang)
plt.title('Diff of Smoothed Phase')
plt.show()
plt.title('Smoothed (ph-phavg)\^2')
plt.plot(errvec)
plt.plot([ires],[errvec[ires]],'ro')
plt.show()
return f0i Example 14
def diff_find_resonance(frec,diffS11,margin=31,doplots=False):
"""
Returns the resonance frequency from maximum of the derivative of the phase
"""
#Smooth data for initial guesses
sReS11 = np.array(smooth(diffS11.real,margin,3))
sImS11 = np.array(smooth(diffS11.imag,margin,3))
sS11 = np.array( [x+1j*y for x,y in zip(sReS11,sImS11) ] )
#Make smoothed phase vector removing 2pi jumps
sArgS11 = np.angle(sS11)
sArgS11 = np.unwrap(sArgS11)
sdiffang = np.diff(sArgS11)
#Get resonance index from maximum of the derivative of the phase
avgph = np.average(sdiffang)
errvec = [np.power(x-avgph,2.) for x in sdiffang]
ires = np.argmax(errvec[margin:-margin])+margin
f0i=frec[ires]
print("Max index: ",ires," Max frequency: ",f0i)
if doplots:
plt.title('Original signal (Re,Im)')
plt.plot(frec,diffS11.real)
plt.plot(frec,diffS11.imag)
plt.plot(frec,np.abs(diffS11))
plt.show()
plt.plot(np.angle(sS11))
plt.title('Smoothed Phase')
plt.axis('auto')
plt.show()
plt.plot(sdiffang[margin:-margin])
plt.title('Diff of Smoothed Phase')
plt.show()
plt.title('Smoothed (ph-phavg)\^2')
plt.plot(errvec)
plt.plot([ires],[errvec[ires]],'ro')
plt.show()
return f0i Example 15
def unwrap_phase_map_0(self):
"""Unwrap order 0 phase map.
Phase is defined modulo pi/2. The Unwrapping is a
reconstruction of the phase so that the distance between two
neighboor pixels is always less than pi/4. Then the real phase
pattern can be recovered and fitted easily.
The idea is the same as with np.unwrap() but in 2D, on a
possibly very noisy map, where a naive 2d unwrapping cannot be
done.
"""
self.phase_map_order_0_unwraped = orb.utils.image.unwrap_phase_map0(
np.copy(self.phase_maps[0]))
# Save unwraped map
phase_map_path = self._get_phase_map_path(0, phase_map_type='unwraped')
self.write_fits(phase_map_path,
orb.cutils.unbin_image(
np.copy(self.phase_map_order_0_unwraped),
self.dimx_unbinned,
self.dimy_unbinned),
fits_header=self._get_phase_map_header(
0, phase_map_type='unwraped'),
overwrite=self.overwrite)
if self.indexer is not None:
self.indexer['phase_map_unwraped_0'] = phase_map_path Example 16
def main():
filename = 'posegraph.posegraph'
plot = False
windowLength = 11
windowType = 'hanning'
if not os.path.isfile(filename):
print 'could not find posegraph file:', filename
sys.exit(1)
data = np.loadtxt(filename)
poseTimes = data[:,0]
poses = np.array(data[:,1:])
poses = np.array([to_xyzrpy(pose) for pose in poses])
poses[:,3:] = np.unwrap(poses[:,3:])
for i in range(poses.shape[1]):
poses[:,i] = smooth(poses[:,i], window_len=windowLength, window=windowType)
if plot:
plot(poses, poseTimes)
poses = np.array([to_xyzquat(pose) for pose in poses])
poses = np.hstack([np.reshape(poseTimes, (-1,1)), poses])
np.savetxt('posegraph_smoothed.posegraph', poses) Example 17
def __init__(self, real_signal, sampling):
self._fs = 1/sampling
self._analytic_signal = hilbert(real_signal)
self._amplitude_envelope = np.abs(self._analytic_signal)
self._instantaneous_phase = np.unwrap(np.angle(self._analytic_signal))
self._instantaneous_frequency = (np.diff(self._instantaneous_phase) /
(2.0*np.pi) * self._fs)
self._instantaneous_frequency = np.insert(self._instantaneous_frequency, 0, np.nan) Example 18
def unwrap(p, discont=pi, axis=-1):
"""
Unwrap by changing deltas between values to 2*pi complement.
Unwrap radian phase `p` by changing absolute jumps greater than
`discont` to their 2*pi complement along the given axis.
Parameters
----------
p : array_like
Input array.
discont : float, optional
Maximum discontinuity between values, default is ``pi``.
axis : int, optional
Axis along which unwrap will operate, default is the last axis.
Returns
-------
out : ndarray
Output array.
See Also
--------
rad2deg, deg2rad
Notes
-----
If the discontinuity in `p` is smaller than ``pi``, but larger than
`discont`, no unwrapping is done because taking the 2*pi complement
would only make the discontinuity larger.
Examples
--------
>>> phase = np.linspace(0, np.pi, num=5)
>>> phase[3:] += np.pi
>>> phase
array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531])
>>> np.unwrap(phase)
array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ])
"""
p = asarray(p)
nd = len(p.shape)
dd = diff(p, axis=axis)
slice1 = [slice(None, None)]*nd # full slices
slice1[axis] = slice(1, None)
ddmod = mod(dd + pi, 2*pi) - pi
_nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
ph_correct = ddmod - dd
_nx.copyto(ph_correct, 0, where=abs(dd) < discont)
up = array(p, copy=True, dtype='d')
up[slice1] = p[slice1] + ph_correct.cumsum(axis)
return up Example 19
def unwrap(p, discont=pi, axis=-1):
"""
Unwrap by changing deltas between values to 2*pi complement.
Unwrap radian phase `p` by changing absolute jumps greater than
`discont` to their 2*pi complement along the given axis.
Parameters
----------
p : array_like
Input array.
discont : float, optional
Maximum discontinuity between values, default is ``pi``.
axis : int, optional
Axis along which unwrap will operate, default is the last axis.
Returns
-------
out : ndarray
Output array.
See Also
--------
rad2deg, deg2rad
Notes
-----
If the discontinuity in `p` is smaller than ``pi``, but larger than
`discont`, no unwrapping is done because taking the 2*pi complement
would only make the discontinuity larger.
Examples
--------
>>> phase = np.linspace(0, np.pi, num=5)
>>> phase[3:] += np.pi
>>> phase
array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531])
>>> np.unwrap(phase)
array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ])
"""
p = asarray(p)
nd = len(p.shape)
dd = diff(p, axis=axis)
slice1 = [slice(None, None)]*nd # full slices
slice1[axis] = slice(1, None)
ddmod = mod(dd + pi, 2*pi) - pi
_nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
ph_correct = ddmod - dd
_nx.copyto(ph_correct, 0, where=abs(dd) < discont)
up = array(p, copy=True, dtype='d')
up[slice1] = p[slice1] + ph_correct.cumsum(axis)
return up Example 20
def analytic_signal(signal, fb, fs=128, order=3):
""" Passband filtering and Hilbert transformation
Parameters
----------
signal: real array-like, shape(n_channels, n_samples)
Input signal
fb: list of length 2
The low and high frequencies
fs: int
Sampling frequency
order : int
Filter order
Returns
-------
filtered_signal: real array-like, shape(n_channels, n_samples)
The input signal, filtered within the given frequencies
hilberted_signal: complex array-like, shape(n_channels, n_samples)
The Hilbert representation of the input signal
unwrapped_phase: real array-like, shape(n_channels, n_samples)
The unwrapped phase of the Hilbert representation
Notes
-----
Internally, we use SciPy's Butterworth implementation (`scipy.signal.butter`)
and the two-pass filter `scipy.signal.filtfilt` to achieve results identical
to MATLAB.
"""
fs = float(fs)
passband = [fb[0] / (fs / 2.0), fb[1] / (fs / 2.0)]
passband = np.ravel(passband)
b, a = scipy.signal.butter(
order, passband, 'bandpass', analog=False, output='ba')
filtered_signal = scipy.signal.filtfilt(b, a, signal)
hilberted_signal = scipy.signal.hilbert(filtered_signal)
unwrapped_phase = np.unwrap(np.angle(hilberted_signal))
return (filtered_signal, hilberted_signal, unwrapped_phase) Example 21
def test_temporal_learning(flen=[5]):
'''
Function that computes CSP filters then uses those with the temporal filter MTL
idea, and confirms that the output has a spectral profile that is similar to expected.
Generate y values from trace of filtered covariance to ensure that is not an issue
'''
def covmat(x,y):
return (1/(x.shape[1]-1)*x.dot(y.T))
D = scio.loadmat('/is/ei/vjayaram/code/python/pyMTL_MunichMI.mat')
data = D['T'].ravel()
labels = D['Y'].ravel()
fig = plt.figure()
fig.suptitle('Recovered frequency filters for various filter lengths')
model = TemporalBRC(max_prior_iter=100)
# compute cross-covariance matrices and generate X
for find,freq in enumerate(flen):
X = []
for tind,d in enumerate(data):
d = d.transpose((2,0,1))
x = np.zeros((d.shape[0], freq))
nsamples = d.shape[2]
for ind, t in enumerate(d):
for j in range(freq):
C = covmat(t[:,0:(nsamples-j)],t[:,j:])
x[ind,j] = np.trace(C + C.T)
X.append(x)
labels[tind] = labels[tind].ravel()
model.fit_multi_task(X,labels)
b = numerics.solve_fir_coef(model.prior.mu.flatten())[0]
# look at filter properties
w,h = freqz(b[1:],worN=100)
w = w*500/2/np.pi # convert to Hz
ax1 = fig.add_subplot(3,3,find+1)
plt.plot(w, 20 * np.log10(abs(h)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [Hz]')
ax2 = ax1.twinx()
angles = np.unwrap(np.angle(h))
plt.plot(w, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.grid()
plt.axis('tight')
plt.show() Example 22
def GetAllData(self,keep_uncal=True):
pars,parnames = self.GetTraceNames()
self.SetActiveTrace(pars[0])
names = ['Frequency (Hz)']
alltrc = [self.GetFrequency()]
for pp in parnames:
names.append('%sre ()' % pp)
names.append('%sim ()' % pp)
names.append('%sdB (dB)' % pp)
names.append('%sPh (rad)' % pp)
for par in pars:
self.SetActiveTrace(par)
yyre,yyim = self.GetTraceData()
alltrc.append(yyre)
alltrc.append(yyim)
yydb = 20.*np.log10( np.abs(yyre+1j*yyim) )
yyph = np.unwrap(np.angle(yyre+1j*yyim))
alltrc.append(yydb)
alltrc.append(yyph)
Cal = self.GetCal()
if Cal and keep_uncal:
for pp in parnames:
names.append('%sre unc ()' % pp)
names.append('%sim unc ()' % pp)
names.append('%sdB unc (dB)' % pp)
names.append('%sPh unc (rad)' % pp)
self.CalOff()
for par in pars:
self.SetActiveTrace(par)
yyre,yyim = self.GetTraceData()
alltrc.append(yyre)
alltrc.append(yyim)
yydb = 20.*np.log10( np.abs(yyre+1j*yyim) )
yyph = np.unwrap(np.angle(yyre+1j*yyim))
alltrc.append(yydb)
alltrc.append(yyph)
self.CalOn()
final = stlabdict()
for name,data in zip(names,alltrc):
final[name]=data
final.addparcolumn('Power (dBm)',self.GetPower())
return final Example 23
def unwrap(p, discont=pi, axis=-1):
"""
Unwrap by changing deltas between values to 2*pi complement.
Unwrap radian phase `p` by changing absolute jumps greater than
`discont` to their 2*pi complement along the given axis.
Parameters
----------
p : array_like
Input array.
discont : float, optional
Maximum discontinuity between values, default is ``pi``.
axis : int, optional
Axis along which unwrap will operate, default is the last axis.
Returns
-------
out : ndarray
Output array.
See Also
--------
rad2deg, deg2rad
Notes
-----
If the discontinuity in `p` is smaller than ``pi``, but larger than
`discont`, no unwrapping is done because taking the 2*pi complement
would only make the discontinuity larger.
Examples
--------
>>> phase = np.linspace(0, np.pi, num=5)
>>> phase[3:] += np.pi
>>> phase
array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531])
>>> np.unwrap(phase)
array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ])
"""
p = asarray(p)
nd = len(p.shape)
dd = diff(p, axis=axis)
slice1 = [slice(None, None)]*nd # full slices
slice1[axis] = slice(1, None)
ddmod = mod(dd + pi, 2*pi) - pi
_nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
ph_correct = ddmod - dd
_nx.copyto(ph_correct, 0, where=abs(dd) < discont)
up = array(p, copy=True, dtype='d')
up[slice1] = p[slice1] + ph_correct.cumsum(axis)
return up Example 24
def unwrap(p, discont=pi, axis=-1):
"""
Unwrap by changing deltas between values to 2*pi complement.
Unwrap radian phase `p` by changing absolute jumps greater than
`discont` to their 2*pi complement along the given axis.
Parameters
----------
p : array_like
Input array.
discont : float, optional
Maximum discontinuity between values, default is ``pi``.
axis : int, optional
Axis along which unwrap will operate, default is the last axis.
Returns
-------
out : ndarray
Output array.
See Also
--------
rad2deg, deg2rad
Notes
-----
If the discontinuity in `p` is smaller than ``pi``, but larger than
`discont`, no unwrapping is done because taking the 2*pi complement
would only make the discontinuity larger.
Examples
--------
>>> phase = np.linspace(0, np.pi, num=5)
>>> phase[3:] += np.pi
>>> phase
array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531])
>>> np.unwrap(phase)
array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ])
"""
p = asarray(p)
nd = len(p.shape)
dd = diff(p, axis=axis)
slice1 = [slice(None, None)]*nd # full slices
slice1[axis] = slice(1, None)
ddmod = mod(dd + pi, 2*pi) - pi
_nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
ph_correct = ddmod - dd
_nx.copyto(ph_correct, 0, where=abs(dd) < discont)
up = array(p, copy=True, dtype='d')
up[slice1] = p[slice1] + ph_correct.cumsum(axis)
return up Example 25
def unwrap(p, discont=pi, axis=-1):
"""
Unwrap by changing deltas between values to 2*pi complement.
Unwrap radian phase `p` by changing absolute jumps greater than
`discont` to their 2*pi complement along the given axis.
Parameters
----------
p : array_like
Input array.
discont : float, optional
Maximum discontinuity between values, default is ``pi``.
axis : int, optional
Axis along which unwrap will operate, default is the last axis.
Returns
-------
out : ndarray
Output array.
See Also
--------
rad2deg, deg2rad
Notes
-----
If the discontinuity in `p` is smaller than ``pi``, but larger than
`discont`, no unwrapping is done because taking the 2*pi complement
would only make the discontinuity larger.
Examples
--------
>>> phase = np.linspace(0, np.pi, num=5)
>>> phase[3:] += np.pi
>>> phase
array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531])
>>> np.unwrap(phase)
array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ])
"""
p = asarray(p)
nd = len(p.shape)
dd = diff(p, axis=axis)
slice1 = [slice(None, None)]*nd # full slices
slice1[axis] = slice(1, None)
ddmod = mod(dd + pi, 2*pi) - pi
_nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
ph_correct = ddmod - dd
_nx.copyto(ph_correct, 0, where=abs(dd) < discont)
up = array(p, copy=True, dtype='d')
up[slice1] = p[slice1] + ph_correct.cumsum(axis)
return up Example 26
def STFT(x, window_size, noverlap, time_start, Ts, mode='psd'):
print 'Ts:', Ts
f0 = 1.0/Ts
print 'F Zero:', f0
print 'Frquencia de Nyquist:', f0/2, 'Hz'
print 'Resolução em frquencia:', f0/window_size, 'Hz'
print 'Resolução temporal:', Ts, 'seconds'
#
window = scipy.hanning(window_size)
stft_data = np.array([np.fft.fft(window*data)
for data in SlidingWindow(x, len(window), noverlap)]
)
#
if len(window) % 2:
numFreqs = stft_data.shape[1]/2+1
else:
numFreqs = stft_data.shape[1]/2
freqs = np.fft.fftfreq(window_size, Ts)[:numFreqs]
stft_data = stft_data[:, :numFreqs]
time_values = np.arange(window_size/2,
len(x)-window_size/2+1,
window_size-noverlap) * Ts + time_start
#
if mode == 'psd':
# Also include scaling factors for one-sided densities and dividing by
# the sampling frequency, if desired. Scale everything, except the DC
# component and the NFFT/2 component:
stft_data *= 2.
# MATLAB divides by the sampling frequency so that density function
# has units of dB/Hz and can be integrated by the plotted frequency
# values. Perform the same scaling here.
#if scale_by_freq:
scale_by_freq = True
if scale_by_freq:
Fs = 1/Ts
stft_data /= Fs
# Scale the spectrum by the norm of the window to compensate for
# windowing loss; see Bendat & Piersol Sec 11.5.2.
stft_data /= (np.abs(window)**2).sum()
else:
# In this case, preserve power in the segment, not amplitude
stft_data /= np.abs(window).sum()**2
stft_data = stft_data * np.conjugate(stft_data)
stft_data = stft_data.real
elif mode == 'magnitude':
stft_data = np.absolute(stft_data)
elif mode == 'angle' or mode == 'phase':
stft_data = np.angle(stft_data)
if mode == 'phase':
stft_data = np.unwrap(stft_data, axis=0)
return stft_data.T, freqs, time_values Example 27
def unwrap(p, discont=pi, axis=-1):
"""
Unwrap by changing deltas between values to 2*pi complement.
Unwrap radian phase `p` by changing absolute jumps greater than
`discont` to their 2*pi complement along the given axis.
Parameters
----------
p : array_like
Input array.
discont : float, optional
Maximum discontinuity between values, default is ``pi``.
axis : int, optional
Axis along which unwrap will operate, default is the last axis.
Returns
-------
out : ndarray
Output array.
See Also
--------
rad2deg, deg2rad
Notes
-----
If the discontinuity in `p` is smaller than ``pi``, but larger than
`discont`, no unwrapping is done because taking the 2*pi complement
would only make the discontinuity larger.
Examples
--------
>>> phase = np.linspace(0, np.pi, num=5)
>>> phase[3:] += np.pi
>>> phase
array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531])
>>> np.unwrap(phase)
array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ])
"""
p = asarray(p)
nd = len(p.shape)
dd = diff(p, axis=axis)
slice1 = [slice(None, None)]*nd # full slices
slice1[axis] = slice(1, None)
ddmod = mod(dd + pi, 2*pi) - pi
_nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
ph_correct = ddmod - dd
_nx.copyto(ph_correct, 0, where=abs(dd) < discont)
up = array(p, copy=True, dtype='d')
up[slice1] = p[slice1] + ph_correct.cumsum(axis)
return up 