def est_tone_phase(sdata,a,f,sr):
samples = len(sdata)
points = 360
rms = numpy.zeros(points)
sum_min = numpy.sum(numpy.square(sdata))
min_index = 0
for offset in xrange(points):
sum = 0
phase = pi*offset/180.0
for i in xrange(samples):
diff = (sdata[i] - a*cos(2*pi*i*f/(sr/2.0) + phase))
sum += diff*diff
rms[offset] = sum
if (sum < sum_min):
sum_min = sum
min_index = offset
#print "sum_min",sum_min,' index = ',min_index
min_phase = pi*(min_index)/180.0
#print "min for phase sweep is ",sum_min,' at offset ',min_index
return min_phase
python类phase()的实例源码
def old_est_tone_phase(sdata,a,f,sr):
samples = len(sdata)
points = 360
rms = numpy.zeros(points)
sum_min = numpy.sum(numpy.square(sdata))
min_index = 0
for offset in xrange(points):
sum = 0
phase = pi*offset/180.0
for i in xrange(samples):
diff = (sdata[i] - a*cos(2*pi*i*f/sr + phase))
sum += diff*diff
rms[offset] = sum
if (sum < sum_min):
sum_min = sum
min_index = offset
#print "sum_min",sum_min,' index = ',min_index
min_phase = pi*(min_index)/180.0
#print "min for phase sweep is ",sum_min,' at offset ',min_index
return min_phase
def find_min_phase(sdata,a,f,sr,phase):
rms1 = 0
rms2 = 0
rms3 = 0
samples = len(sdata)
for i in xrange(samples):
diff1 = (sdata[i] - a*cos(2*pi*i*f/sr + phase[0]))
rms1 += diff1*diff1
diff2 = (sdata[i] - a*cos(2*pi*i*f/sr + phase[1]))
rms2 += diff2*diff2
diff3 = (sdata[i] - a*cos(2*pi*i*f/sr + phase[2]))
rms3 += diff3*diff3
rms = numpy.zeros(3)
rms[0] = rms1
rms[1] = rms2
rms[2] = rms3
i = numpy.argmin(rms)
p = phase[i]
return i,p
def est_tone_phase(sdata,a,f,sr):
delta = 120
min_ang = 0.5
p = 0
phase = numpy.zeros(3)
phase[0] = p+(-delta/180.0)*pi
phase[1] = p
phase[2] = p+(delta/180.0)*pi
while (delta > min_ang):
(i,p) = find_min_phase(sdata,a,f,sr,phase)
delta = delta/2.0
phase[0] = p+(-delta/180.0)*pi
phase[1] = p
phase[2] = p+(delta/180.0)*pi
#print "p = ",(180.0*p/pi),'delta = ',delta
min_phase = p
#print "min for phase sweep is ",sum_min,' at offset ',min_index
return min_phase
def arc_to(self, x1, y1):
x0, y0 = self.get_current_point()
dx, dy = x1 - x0, y1 - y0
if abs(dx) < 1e-8:
self.line_to(x1, y1)
else:
center = 0.5 * (x0 + x1) + 0.5 * (y0 + y1) * dy / dx
theta0 = cmath.phase(x0 - center + y0*1j)
theta1 = cmath.phase(x1 - center + y1*1j)
r = abs(x0 - center + y0*1j)
# we must ensure that the arc ends at (x1, y1)
if x0 < x1:
self.arc_negative(center, 0, r, theta0, theta1)
else:
self.arc(center, 0, r, theta0, theta1)
def phase2t(self, psi):
"""Given phase -pi < psi <= pi,
returns the t value such that
exp(1j*psi) = self.u1transform(self.point(t)).
"""
def _deg(rads, domain_lower_limit):
# Convert rads to degrees in [0, 360) domain
degs = degrees(rads % (2*pi))
# Convert to [domain_lower_limit, domain_lower_limit + 360) domain
k = domain_lower_limit // 360
degs += k * 360
if degs < domain_lower_limit:
degs += 360
return degs
if self.delta > 0:
degs = _deg(psi, domain_lower_limit=self.theta)
else:
degs = _deg(psi, domain_lower_limit=self.theta)
return (degs - self.theta)/self.delta
def tone_est(sdata,sr):
samples = len(sdata)
fft_size = 2**int(floor(log(samples)/log(2.0)))
freq = fft(sdata[0:fft_size])
pdata = numpy.zeros(fft_size)
for i in xrange(fft_size): pdata[i] = abs(freq[i])
peak = 0
peak_index = 0
for i in xrange(fft_size/2):
if (pdata[i] > peak):
peak = pdata[i]
peak_index = i
R = peak*peak;
p = (freq[peak_index+1].real * freq[peak_index].real + freq[peak_index+1].imag * freq[peak_index].imag)/R
g = -p/(1.0-p)
q = (freq[peak_index-1].real * freq[peak_index].real + freq[peak_index-1].imag * freq[peak_index].imag)/R
e = q/(1.0-q)
if ((p>0) and (q>0)):
d = p
else:
d = q
u = peak_index + d
print "peak is at ",peak_index,"(",u,") and is ",peak
#print "u = ",0.5*u*sr/fft_size,' f[0] = ',f[0]
sum_phase = freq[peak_index-1]*c(-1,d) + freq[peak_index]*c(0,d) + freq[peak_index+1]*c(1,d)
sum_c_sq = abs(c(-1,d))*abs(c(-1,d)) + abs(c(0,d))*abs(c(0,d)) + abs(c(1,d))*abs(c(1,d))
amp = (abs(sum_phase)/sum_c_sq)/fft_size
phase_r = cmath.phase(sum_phase)
freq_est = 0.5*u*sr/fft_size
return (amp,freq_est,phase_r)
def tone_est_near_index(sdata,index,range,sr):
samples = len(sdata)
fft_size = 2**int(floor(log(samples)/log(2.0)))
freq = fft(sdata[0:fft_size])
pdata = numpy.zeros(fft_size)
for i in xrange(fft_size): pdata[i] = abs(freq[i])
peak = 0
peak_index = 0
for i in xrange(2*range):
if (pdata[index+i-range] > peak):
peak = pdata[index+i-range]
peak_index = index+i-range
print "peak is at ",peak_index," and is ",peak
R = peak*peak;
p = (freq[peak_index+1].real * freq[peak_index].real + freq[peak_index+1].imag * freq[peak_index].imag)/R
g = -p/(1.0-p)
q = (freq[peak_index-1].real * freq[peak_index].real + freq[peak_index-1].imag * freq[peak_index].imag)/R
e = q/(1.0-q)
if ((p>0) and (q>0)):
d = p
else:
d = q
u = peak_index + d
sum_phase = freq[peak_index-1]*c(-1,d) + freq[peak_index]*c(0,d) + freq[peak_index+1]*c(1,d)
sum_c_sq = abs(c(-1,d))*abs(c(-1,d)) + abs(c(0,d))*abs(c(0,d)) + abs(c(1,d))*abs(c(1,d))
amp = (abs(sum_phase)/sum_c_sq)/fft_size
phase_r = cmath.phase(sum_phase)
freq_est = 0.5*u*sr/fft_size
return (amp,freq_est,phase_r)
def tone_est(sdata,sr):
samples = len(sdata)
fft_size = 2**int(floor(log(samples)/log(2.0)))
freq = fft(sdata[0:fft_size])
pdata = numpy.zeros(fft_size)
for i in xrange(fft_size): pdata[i] = abs(freq[i])
peak = 0
peak_index = 0
for i in xrange(fft_size/2):
if (pdata[i] > peak):
peak = pdata[i]
peak_index = i
R = peak*peak;
p = (freq[peak_index+1].real * freq[peak_index].real + freq[peak_index+1].imag * freq[peak_index].imag)/R
q = (freq[peak_index-1].real * freq[peak_index].real + freq[peak_index-1].imag * freq[peak_index].imag)/R
g = -p/(1.0-p)
e = q/(1.0-q)
if ((p>0) and (q>0)):
d = p
else:
d = q
u = peak_index + d
freq_est = u*sr/fft_size
if (debug_estimates):
print "peak is at ",peak_index,"(",u,") and is ",peak
#d = 0.5*(p+q) + h(p*p) + h(q*q)
#print "other peak index (2)", u+d
sum_phase = freq[peak_index-1]*c(-1,d) + freq[peak_index]*c(0,d) + freq[peak_index+1]*c(1,d)
sum_c_sq = abs(c(-1,d))*abs(c(-1,d)) + abs(c(0,d))*abs(c(0,d)) + abs(c(1,d))*abs(c(1,d))
amp = (abs(sum_phase)/sum_c_sq)/fft_size
phase_r = cmath.phase(sum_phase)
return (amp,freq_est,phase_r)
def tone_est_above_index(sdata,index,sr):
samples = len(sdata)
fft_size = 2**int(floor(log(samples)/log(2.0)))
freq = fft(sdata[0:fft_size])
pdata = numpy.zeros(fft_size)
for i in xrange(fft_size): pdata[i] = abs(freq[i])
peak = 0
peak_index = 0
for i in xrange(fft_size/2):
if (i > index):
if (pdata[i] > peak):
peak = pdata[i]
peak_index = i
R = peak*peak;
p = (freq[peak_index+1].real * freq[peak_index].real + freq[peak_index+1].imag * freq[peak_index].imag)/R
g = -p/(1.0-p)
q = (freq[peak_index-1].real * freq[peak_index].real + freq[peak_index-1].imag * freq[peak_index].imag)/R
e = q/(1.0-q)
if ((p>0) and (q>0)):
d = p
else:
d = q
u = peak_index + d
if (debug_estimates):
print "peak is at ",peak_index,"(",u,") and is ",peak
sum_phase = freq[peak_index-1]*c(-1,d) + freq[peak_index]*c(0,d) + freq[peak_index+1]*c(1,d)
sum_c_sq = abs(c(-1,d))*abs(c(-1,d)) + abs(c(0,d))*abs(c(0,d)) + abs(c(1,d))*abs(c(1,d))
amp = (abs(sum_phase)/sum_c_sq)/fft_size
phase_r = cmath.phase(sum_phase)
freq_est = u*sr/fft_size
return (amp,freq_est,phase_r)
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assertTrue(math.isnan(phase(z)))
def mean_heading(headings):
"""
Calculate the average heading from a list of headings
:param headings:
:return: average heading
"""
vectors = [cmath.rect(1, radians(angle)) for angle in headings]
vector_sum = sum(vectors)
return degrees(cmath.phase(vector_sum))
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assertTrue(math.isnan(phase(z)))
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assertTrue(math.isnan(phase(z)))
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assertTrue(math.isnan(phase(z)))
def mean_angle(self, deg):
return (phase(sum(rect(1, d) for d in deg)/len(deg)))
def fm_detect(X, prev, shift):
vals = array.array('f')
for x in X:
y = shift + cmath.phase(x * prev.conjugate()) / math.pi
vals.append(y)
prev = x
return vals
def process(self, samples):
Y = array.array('f')
prev = self._prev
for x in samples:
y = cmath.phase(x * prev.conjugate()) / math.pi
Y.append(y)
prev = x
self._prev = prev
return Y
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assertTrue(math.isnan(phase(z)))
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assertTrue(math.isnan(phase(z)))
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assertTrue(math.isnan(phase(z)))
def parse_DFT_output(FmList, threshold=0.001):
outputs = []
for (i, Fm) in enumerate(FmList):
if abs(Fm) > threshold:
frequency = i
amplitude = abs(Fm) * 2.0
phase_angle = int(((cmath.phase(Fm) + pi2 + pi2 / 4.0) % pi2) / pi2 * 360 + 0.5)
outputs.append((amplitude, frequency, phase_angle))
return outputs
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assertTrue(math.isnan(phase(z)))
def curvature(self, t):
"""returns the curvature of the segment at t."""
return segment_curvature(self, t)
# def icurvature(self, kappa):
# """returns a list of t-values such that 0 <= t<= 1 and
# seg.curvature(t) = kappa."""
#
# a, b = self.radius.real, self.radius.imag
# if kappa > min(a, b)/max(a, b)**2 or kappa <= 0:
# return []
# if a==b:
# if kappa != 1/a:
# return []
# else:
# raise ValueError(
# "The .icurvature() method for Arc elements with "
# "radius.real == radius.imag (i.e. circle segments) "
# "will raise this exception when kappa is 1/radius.real as "
# "this is true at every point on the circle segment.")
#
# # kappa = a*b / (a^2sin^2(tau) + b^2cos^2(tau))^(3/2), tau=2*pi*phase
# sin2 = np.poly1d([1, 0])
# p = kappa**2*(a*sin2 + b*(1 - sin2))**3 - a*b
# sin2s = polyroots01(p)
# taus = []
#
# for sin2 in sin2s:
# taus += [np.arcsin(sqrt(sin2)), np.arcsin(-sqrt(sin2))]
#
# # account for the other branch of arcsin
# sgn = lambda x: x/abs(x) if x else 0
# other_taus = [sgn(tau)*np.pi - tau for tau in taus if abs(tau) != np.pi/2]
# taus = taus + other_taus
#
# # get rid of points not included in segment
# ts = [phase2t(tau) for tau in taus]
#
# return [t for t in ts if 0<=t<=1]
def curvature(self, t):
"""returns the curvature of the segment at t."""
return segment_curvature(self, t)
# def icurvature(self, kappa):
# """returns a list of t-values such that 0 <= t<= 1 and
# seg.curvature(t) = kappa."""
#
# a, b = self.radius.real, self.radius.imag
# if kappa > min(a, b)/max(a, b)**2 or kappa <= 0:
# return []
# if a==b:
# if kappa != 1/a:
# return []
# else:
# raise ValueError(
# "The .icurvature() method for Arc elements with "
# "radius.real == radius.imag (i.e. circle segments) "
# "will raise this exception when kappa is 1/radius.real as "
# "this is true at every point on the circle segment.")
#
# # kappa = a*b / (a^2sin^2(tau) + b^2cos^2(tau))^(3/2), tau=2*pi*phase
# sin2 = np.poly1d([1, 0])
# p = kappa**2*(a*sin2 + b*(1 - sin2))**3 - a*b
# sin2s = polyroots01(p)
# taus = []
#
# for sin2 in sin2s:
# taus += [np.arcsin(sqrt(sin2)), np.arcsin(-sqrt(sin2))]
#
# # account for the other branch of arcsin
# sgn = lambda x: x/abs(x) if x else 0
# other_taus = [sgn(tau)*np.pi - tau for tau in taus if abs(tau) != np.pi/2]
# taus = taus + other_taus
#
# # get rid of points not included in segment
# ts = [phase2t(tau) for tau in taus]
#
# return [t for t in ts if 0<=t<=1]
def tone_est_above_index(sdata,index,sr):
samples = len(sdata)
fft_size = 2**int(floor(log(samples)/log(2.0)))
freq = fft(sdata[0:fft_size])
pdata = numpy.zeros(fft_size)
for i in xrange(fft_size): pdata[i] = abs(freq[i])
peak = 0
peak_index = 0
for i in xrange(fft_size/2):
if (i > index):
if (pdata[i] > peak):
peak = pdata[i]
peak_index = i
print "peak is at ",peak_index," and is ",peak
R = peak*peak;
p = (freq[peak_index+1].real * freq[peak_index].real + freq[peak_index+1].imag * freq[peak_index].imag)/R
g = -p/(1.0-p)
q = (freq[peak_index-1].real * freq[peak_index].real + freq[peak_index-1].imag * freq[peak_index].imag)/R
e = q/(1.0-q)
if ((p>0) and (q>0)):
d = p
else:
d = q
u = peak_index + d
sum_phase = freq[peak_index-1]*c(-1,d) + freq[peak_index]*c(0,d) + freq[peak_index+1]*c(1,d)
sum_c_sq = abs(c(-1,d))*abs(c(-1,d)) + abs(c(0,d))*abs(c(0,d)) + abs(c(1,d))*abs(c(1,d))
amp = (abs(sum_phase)/sum_c_sq)/fft_size
phase_r = cmath.phase(sum_phase)
freq_est = 0.5*u*sr/fft_size
return (amp,freq_est,phase_r)
# REMOVAL CASES
def intersect(self, other_seg, tol=1e-12):
"""NOT FULLY IMPLEMENTED. Finds the intersections of two segments.
returns a list of tuples (t1, t2) such that
self.point(t1) == other_seg.point(t2).
Note: This will fail if the two segments coincide for more than a
finite collection of points.
Note: Arc related intersections are only partially supported, i.e. are
only half-heartedly implemented and not well tested. Please feel free
to let me know if you're interested in such a feature -- or even better
please submit an implementation if you want to code one."""
if is_bezier_segment(other_seg):
u1poly = self.u1transform(other_seg.poly())
u1poly_mag2 = real(u1poly)**2 + imag(u1poly)**2
t2s = polyroots01(u1poly_mag2 - 1)
t1s = [self.phase2t(phase(u1poly(t2))) for t2 in t2s]
return list(zip(t1s, t2s))
elif isinstance(other_seg, Arc):
assert other_seg != self
# This could be made explicit to increase efficiency
longer_length = max(self.length(), other_seg.length())
inters = bezier_intersections(self, other_seg,
longer_length=longer_length,
tol=tol, tol_deC=tol)
# ad hoc fix for redundant solutions
if len(inters) > 2:
def keyfcn(tpair):
t1, t2 = tpair
return abs(self.point(t1) - other_seg.point(t2))
inters.sort(key=keyfcn)
for idx in range(1, len(inters)-1):
if (abs(inters[idx][0] - inters[idx + 1][0])
< abs(inters[idx][0] - inters[0][0])):
return [inters[0], inters[idx]]
else:
return [inters[0], inters[-1]]
return inters
else:
raise TypeError("other_seg should be a Arc, Line, "
"QuadraticBezier, or CubicBezier object.")
def radialrange(self, origin, return_all_global_extrema=False):
"""returns the tuples (d_min, t_min) and (d_max, t_max) which minimize
and maximize, respectively, the distance,
d = |self.point(t)-origin|."""
u1orig = self.u1transform(origin)
if abs(u1orig) == 1: # origin lies on ellipse
t = self.phase2t(phase(u1orig))
d_min = 0
# Transform to a coordinate system where the ellipse is centered
# at the origin and its axes are horizontal/vertical
zeta0 = self.centeriso(origin)
a, b = self.radius.real, self.radius.imag
x0, y0 = zeta0.real, zeta0.imag
# Find t s.t. z'(t)
a2mb2 = (a**2 - b**2)
if u1orig.imag: # x != x0
coeffs = [a2mb2**2,
2*a2mb2*b**2*y0,
(-a**4 + (2*a**2 - b**2 + y0**2)*b**2 + x0**2)*b**2,
-2*a2mb2*b**4*y0,
-b**6*y0**2]
ys = polyroots(coeffs, realroots=True,
condition=lambda r: -b <= r <= b)
xs = (a*sqrt(1 - y**2/b**2) for y in ys)
ts = [self.phase2t(phase(self.u1transform(self.icenteriso(
complex(x, y))))) for x, y in zip(xs, ys)]
else: # This case is very similar, see notes and assume instead y0!=y
b2ma2 = (b**2 - a**2)
coeffs = [b2ma2**2,
2*b2ma2*a**2*x0,
(-b**4 + (2*b**2 - a**2 + x0**2)*a**2 + y0**2)*a**2,
-2*b2ma2*a**4*x0,
-a**6*x0**2]
xs = polyroots(coeffs, realroots=True,
condition=lambda r: -a <= r <= a)
ys = (b*sqrt(1 - x**2/a**2) for x in xs)
ts = [self.phase2t(phase(self.u1transform(self.icenteriso(
complex(x, y))))) for x, y in zip(xs, ys)]
raise _NotImplemented4ArcException
def intersect(self, other_seg, tol=1e-12):
"""NOT FULLY IMPLEMENTED. Finds the intersections of two segments.
returns a list of tuples (t1, t2) such that
self.point(t1) == other_seg.point(t2).
Note: This will fail if the two segments coincide for more than a
finite collection of points.
Note: Arc related intersections are only partially supported, i.e. are
only half-heartedly implemented and not well tested. Please feel free
to let me know if you're interested in such a feature -- or even better
please submit an implementation if you want to code one."""
if is_bezier_segment(other_seg):
u1poly = self.u1transform(other_seg.poly())
u1poly_mag2 = real(u1poly)**2 + imag(u1poly)**2
t2s = polyroots01(u1poly_mag2 - 1)
t1s = [self.phase2t(phase(u1poly(t2))) for t2 in t2s]
return list(zip(t1s, t2s))
elif isinstance(other_seg, Arc):
assert other_seg != self
# This could be made explicit to increase efficiency
longer_length = max(self.length(), other_seg.length())
inters = bezier_intersections(self, other_seg,
longer_length=longer_length,
tol=tol, tol_deC=tol)
# ad hoc fix for redundant solutions
if len(inters) > 2:
def keyfcn(tpair):
t1, t2 = tpair
return abs(self.point(t1) - other_seg.point(t2))
inters.sort(key=keyfcn)
for idx in range(1, len(inters)-1):
if (abs(inters[idx][0] - inters[idx + 1][0])
< abs(inters[idx][0] - inters[0][0])):
return [inters[0], inters[idx]]
else:
return [inters[0], inters[-1]]
return inters
else:
raise TypeError("other_seg should be a Arc, Line, "
"QuadraticBezier, or CubicBezier object.")
def radialrange(self, origin, return_all_global_extrema=False):
"""returns the tuples (d_min, t_min) and (d_max, t_max) which minimize
and maximize, respectively, the distance,
d = |self.point(t)-origin|."""
u1orig = self.u1transform(origin)
if abs(u1orig) == 1: # origin lies on ellipse
t = self.phase2t(phase(u1orig))
d_min = 0
# Transform to a coordinate system where the ellipse is centered
# at the origin and its axes are horizontal/vertical
zeta0 = self.centeriso(origin)
a, b = self.radius.real, self.radius.imag
x0, y0 = zeta0.real, zeta0.imag
# Find t s.t. z'(t)
a2mb2 = (a**2 - b**2)
if u1orig.imag: # x != x0
coeffs = [a2mb2**2,
2*a2mb2*b**2*y0,
(-a**4 + (2*a**2 - b**2 + y0**2)*b**2 + x0**2)*b**2,
-2*a2mb2*b**4*y0,
-b**6*y0**2]
ys = polyroots(coeffs, realroots=True,
condition=lambda r: -b <= r <= b)
xs = (a*sqrt(1 - y**2/b**2) for y in ys)
ts = [self.phase2t(phase(self.u1transform(self.icenteriso(
complex(x, y))))) for x, y in zip(xs, ys)]
else: # This case is very similar, see notes and assume instead y0!=y
b2ma2 = (b**2 - a**2)
coeffs = [b2ma2**2,
2*b2ma2*a**2*x0,
(-b**4 + (2*b**2 - a**2 + x0**2)*a**2 + y0**2)*a**2,
-2*b2ma2*a**4*x0,
-a**6*x0**2]
xs = polyroots(coeffs, realroots=True,
condition=lambda r: -a <= r <= a)
ys = (b*sqrt(1 - x**2/a**2) for x in xs)
ts = [self.phase2t(phase(self.u1transform(self.icenteriso(
complex(x, y))))) for x, y in zip(xs, ys)]
raise _NotImplemented4ArcException
