Python numpy 模块，unwrap() 实例源码
我们从Python开源项目中，提取了以下27个代码示例，用于说明如何使用numpy.unwrap()。
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
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)
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')
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)
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))) ]
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 = {}
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
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
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
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
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)
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)
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
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
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
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)
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)
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
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
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)
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()
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
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
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
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
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
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