Added snippet tests

This commit is contained in:
PhilJ 2022-01-13 23:11:55 -08:00
parent 45fde86060
commit 612d07218e
26 changed files with 237 additions and 7106 deletions

View file

@ -45,3 +45,4 @@ make -C micropython/ports/unix -j${NPROC} axtls
make -C micropython/ports/unix -j${NPROC} USER_C_MODULES="${HERE}" DEBUG=1 STRIP=: MICROPY_PY_FFI=0 MICROPY_PY_BTREE=0 CFLAGS_EXTRA=-DULAB_MAX_DIMS=$dims BUILD=build-$dims PROG=micropython-$dims
bash test-common.sh "${dims}" "micropython/ports/unix/micropython-$dims"

View file

@ -51,6 +51,8 @@ def size(a, axis=None):
return asarray(a).shape[axis]
def nonzero(a):
if not isinstance(a,(np.ndarray)):
a = asarray(a)
x = a.shape
row = x[0]
if len(x) == 1:
@ -65,11 +67,12 @@ def nonzero(a):
for i in range(0,row):
if a[i] != 0:
nonzero_row = numpy.append(nonzero_row,i)
return (nonzero_row,)
return (np.array(nonzero_row, dtype=np.int8),)
for i in range(0,row):
for j in range(0,column):
if a[i,j] != 0:
nonzero_row = numpy.append(nonzero_row,i)
nonzero_col = numpy.append(nonzero_col,j)
return (nonzero_row,nonzero_col)
return (np.array(nonzero_row, dtype=np.int8), np.array(nonzero_col, dtype=np.int8))

View file

@ -6,6 +6,8 @@ def asarray(a, dtype=None):
try:
if dtype is not None:
a = np.array([a], dtype=dtype)
elif isinstance(a, list) or isinstance(a, tuple):
a = np.array(a)
else:
a = np.array([a])
return a

View file

@ -56,3 +56,4 @@ def zeros_like(a, dtype=None, order='K', subok=True, shape=None):
res = np.full(a.shape, 0, dtype=a.dtype)
return res

View file

@ -7,9 +7,9 @@ from ulab import numpy as np
def poly(seq_of_zeros):
seq_of_zeros = atleast_1d(seq_of_zeros)
sh = seq_of_zeros.shape
if len(sh) == 2 and sh[0] == sh[1] and sh[0] != 0:
seq_of_zeros = eigvals(seq_of_zeros)
seq_of_zeros = np.linalg.eig(seq_of_zeros)
elif len(sh) == 1:
dt = seq_of_zeros.dtype
# Let object arrays slip through, e.g. for arbitrary precision

View file

@ -31,12 +31,12 @@ def butter(N, Wn, btype='low', analog=False, output='ba', fs=None):
half-cycles / sample.)
For analog filters, `Wn` is an angular frequency (e.g. rad/s).
btype : {'lowpass', 'highpass', 'bandpass', 'bandstop'}, optional
btype : {'lowpass', 'highpass', 'bandpass', 'bandstop'}, optional
The type of filter. Default is 'lowpass'.
analog : bool, optional
analog : bool, optional
When True, return an analog filter, otherwise a digital filter is
returned.
output : {'ba', 'zpk', 'sos'}, optional
output : {'ba', 'zpk', 'sos'}, optional
Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or
second-order sections ('sos'). Default is 'ba' for backwards
compatibility, but 'sos' should be used for general-purpose filtering.
@ -119,7 +119,6 @@ def butter(N, Wn, btype='low', analog=False, output='ba', fs=None):
return iirfilter(N, Wn, btype=btype, analog=analog,
output=output, ftype='butter', fs=fs)
def iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=False,
ftype='butter', output='ba', fs=None):
"""
@ -305,9 +304,9 @@ def iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=False,
raise ValueError('Must specify a single critical frequency Wn for lowpass or highpass filter')
if btype == 'lowpass':
z, p, k = lp2lp_zpk(z, p, k, wo=warped)
z, p, k = lp2lp_zpk(z, p, k, wo=warped[0])
elif btype == 'highpass':
z, p, k = lp2hp_zpk(z, p, k, wo=warped)
z, p, k = lp2hp_zpk(z, p, k, wo=warped[0])
elif btype in ('bandpass', 'bandstop'):
try:
bw = warped[1] - warped[0]
@ -406,13 +405,11 @@ def _to_complex(a):
result = np.concatenate((result, t), axis=0)
return result
def lexsort(z):
z = _to_tuple(z)
return sorted(range(len(z)), key=lambda i: (z[i][0],abs(z[i][1])))
def do_lexsort(z):
# z = z[np.lexsort((abs(z.imag), z.real))]
def _lexsort(z):
z = _to_tuple(z)
z = sorted(z,key=lambda x:(x[0], abs(x[1])))
return _to_complex(z)
@ -472,11 +469,7 @@ def _cplxreal(z, tol=None):
tol = 100 * abs(7./3 - 4./3 - 1) #np.finfo((1.0 * z).dtype).eps
# Sort by real part, magnitude of imaginary part (speed up further sorting)
#z = z[np.lexsort((abs(z.imag), z.real))]
#z = _to_tuple(z)
#z = sorted(z,key=lambda x:(x[0], abs(x[1])))
#z = _to_complex(z)
z = do_lexsort(z)
z = _lexsort(z)
# Split reals from conjugate pairs
real_indices = abs(z.imag) <= tol * abs(z)
zr = z[real_indices].real
@ -698,28 +691,42 @@ def zpk2sos(z, p, k, pairing='nearest'):
raise ValueError('pairing must be one of %s, not %s'
% (valid_pairings, pairing))
if len(z) == len(p) == 0:
return array([[k, 0., 0., 1., 0., 0.]])
return np.array([[k, 0., 0., 1., 0., 0.]])
# ensure we have the same number of poles and zeros, and make copies
if len(z) != len(p):
p = np.concatenate((p, np.zeros(max(len(z) - len(p), 0), dtype=np.complex)))
z = np.concatenate((z, np.zeros(max(len(p) - len(z), 0), dtype=np.complex)))
if max(len(z) - len(p),0) > 0:
p = np.concatenate((p, np.zeros((max(len(z) - len(p), 0)),dtype=np.complex)))
if max(len(p) - len(z),0) > 0:
z = np.concatenate((z, np.zeros((max(len(p) - len(z), 0)), dtype=np.complex)))
n_sections = (max(len(p), len(z)) + 1) // 2
sos = np.zeros((n_sections, 6))
if len(p) % 2 == 1 and pairing == 'nearest':
p = np.concatenate((p, [0.]))
z = np.concatenate((z, [0.]))
p = np.concatenate((p, np.array([0.],dtype=np.complex)))
z = np.concatenate((z, np.array([0.],dtype=np.complex)))
assert len(p) == len(z)
z = np.array(z, dtype=np.complex)
# Ensure we have complex conjugate pairs
# (note that _cplxreal only gives us one element of each complex pair):
cplx, real = _cplxreal(p)
p = cplx
cplx, real = _cplxreal(z)
z = real
cplx, rel = _cplxreal(p)
if len(rel) > 0 and len(cplx) > 0:
p = np.concatenate((cplx, np.array(rel, dtype=np.complex)))
else:
p = cplx
cplx, rel = _cplxreal(z)
if len(rel) > 0 and len(cplx) > 0:
z = np.concatenate((cplx, np.array(rel, dtype=np.complex)))
else:
z = rel
p_sos = np.zeros((n_sections, 2))
z_sos = zeros_like(p_sos)
@ -730,7 +737,7 @@ def zpk2sos(z, p, k, pairing='nearest'):
p1 = p[p1_idx]
p = delete(p, p1_idx)
# Pair that pole with a zero
if isreal(p1) and np.sum(isreal(p)) == 0:
if isreal(p1) and np.sum([isreal(p)]) == 0:
# Special case to set a first-order section
z1_idx = _nearest_real_complex_idx(z, p1, 'real')
z1 = z[z1_idx]
@ -766,9 +773,12 @@ def zpk2sos(z, p, k, pairing='nearest'):
assert isreal(p2)
else: # real pole, real zero
# pick the next "worst" pole to use
idx = np.nonzero(np.isreal(p))[0]
idx = nonzero(isreal(p))[0]
assert len(idx) > 0
p2_idx = idx[np.argmin(np.abs(abs(p[idx]) - 1))]
a = abs(abs(p[idx[0]]) - 1)
a = np.array([a])
p2_idx = idx[np.argmin(a) - 1]
p2 = p[p2_idx]
# find a real zero to match the added pole
assert isreal(p2)
@ -797,7 +807,6 @@ def zpk2sos(z, p, k, pairing='nearest'):
sos[si] = np.concatenate(x)
return sos
def lp2bp_zpk(z, p, k, wo=1.0, bw=1.0):
r"""
Transform a lowpass filter prototype to a bandpass filter.
@ -855,7 +864,6 @@ def lp2bp_zpk(z, p, k, wo=1.0, bw=1.0):
degree = _relative_degree(z, p)
# Scale poles and zeros to desired bandwidth
# z_lp = z * bw/2
z_lp = []
for x in z:
z_lp.append(x * bw/2)
@ -869,48 +877,23 @@ def lp2bp_zpk(z, p, k, wo=1.0, bw=1.0):
# Square root needs to produce complex result, not NaN
# z_lp = z_lp.astype(complex)
# p_lp = p_lp.astype(complex)
# Duplicate poles and zeros and shift from baseband to +wo and -wo
z_bp = []
for x in z_lp:
z_bp.append(x + np.sqrt(x**2 - wo**2))
for x in z_lp:
z_bp.append(x - np.sqrt(x**2 - wo**2))
z_bp = np.array(z_bp, dtype=np.complex)
# z_bp = np.concatenate((z_lp + np.sqrt(z_lp**2 - wo**2),
# z_lp - np.sqrt(z_lp**2 - wo**2)))
p_bp = []
for x in p_lp:
p_bp.append(x + np.sqrt(x**2 - wo**2))
for x in p_lp:
p_bp.append(x - np.sqrt(x**2 - wo**2))
p_bp = np.array(p_bp, dtype=np.complex)
# p_bp = np.concatenate((p_lp + np.sqrt(p_lp**2 - wo**2),
# p_lp - np.sqrt(p_lp**2 - wo**2)))
if len(z_lp) > 0:
z_bp = np.concatenate((z_lp + np.sqrt(z_lp*z_lp - wo*wo, dtype=np.complex),
z_lp - np.sqrt(z_lp*z_lp - wo*wo, dtype=np.complex)))
z_bp = append(z_bp, np.zeros(degree),dtype=np.complex)
else:
z_bp = np.zeros(degree, dtype=np.complex)
p_bp = np.concatenate((p_lp + np.sqrt(p_lp*p_lp - wo*wo, dtype=np.complex),
p_lp - np.sqrt(p_lp*p_lp - wo*wo, dtype=np.complex)))
# Move degree zeros to origin, leaving degree zeros at infinity for BPF
t = []
for x in z_bp:
t.append(x)
for x in np.zeros(degree):
t.append(x)
z_bp = np.array(t, dtype=np.complex)
# Cancel out gain change from frequency scaling
k_bp = k * bw**degree
@ -965,8 +948,8 @@ def lp2bs_zpk(z, p, k, wo=1.0, bw=1.0):
.. versionadded:: 1.1.0
"""
z = np.atleast_1d(z)
p = np.atleast_1d(p)
z = atleast_1d(z)
p = atleast_1d(p)
wo = float(wo)
bw = float(bw)
@ -976,22 +959,29 @@ def lp2bs_zpk(z, p, k, wo=1.0, bw=1.0):
z_hp = (bw/2) / z
p_hp = (bw/2) / p
# Square root needs to produce complex result, not NaN
z_hp = z_hp.astype(complex)
p_hp = p_hp.astype(complex)
if z_hp == float('inf'):
z_hp = []
# Duplicate poles and zeros and shift from baseband to +wo and -wo
z_bs = concatenate((z_hp + sqrt(z_hp**2 - wo**2),
z_hp - sqrt(z_hp**2 - wo**2)))
p_bs = concatenate((p_hp + sqrt(p_hp**2 - wo**2),
p_hp - sqrt(p_hp**2 - wo**2)))
# Square root needs to produce complex result, not NaN
z_hp = np.array(z_hp, dtype=np.complex)
p_hp = np.array(p_hp, dtype=np.complex)
if len(z_hp) > 0:
# Duplicate poles and zeros and shift from baseband to +wo and -wo
z_bs = np.concatenate((z_hp + np.sqrt(z_hp*z_hp - wo*wo, dtype=np.complex),
z_hp - np.sqrt(z_hp*z_hp - wo*wo, dtype=np.complex)))
else:
z_bs = np.array([], dtype=np.complex)
p_bs = np.concatenate((p_hp + np.sqrt(p_hp*p_hp - wo*wo, dtype=np.complex),
p_hp - np.sqrt(p_hp*p_hp - wo*wo, dtype=np.complex)))
# Move any zeros that were at infinity to the center of the stopband
z_bs = append(z_bs, full(degree, +1j*wo))
z_bs = append(z_bs, full(degree, -1j*wo))
z_bs = append(z_bs, np.full(degree, 1j*wo, dtype=np.complex))
z_bs = append(z_bs, np.full(degree, -1j*wo, dtype=np.complex))
# Cancel out gain change caused by inversion
k_bs = k * real(prod(-z) / prod(-p))
k_bs = k * (prod(-z) / prod(-p)).real
return z_bs, p_bs, k_bs
@ -1056,6 +1046,7 @@ def bilinear_zpk(z, p, k, fs):
>>> plt.ylabel('Magnitude [dB]')
>>> plt.grid()
"""
z = atleast_1d(z)
p = atleast_1d(p)
@ -1070,7 +1061,8 @@ def bilinear_zpk(z, p, k, fs):
# Any zeros that were at infinity get moved to the Nyquist frequency
a = -np.ones(degree) + 0j
z_z = append(z_z, a)
if len(a) > 0:
z_z = append(z_z, a)
# Compensate for gain change
k_z = k * (prod(fs2 - z) / prod(fs2 - p)).real
@ -1215,8 +1207,9 @@ def lp2hp_zpk(z, p, k, wo=1.0):
.. versionadded:: 1.1.0
"""
z = np.atleast_1d(z)
p = np.atleast_1d(p)
z = atleast_1d(z)
p = atleast_1d(p)
wo = float(wo)
degree = _relative_degree(z, p)
@ -1226,11 +1219,14 @@ def lp2hp_zpk(z, p, k, wo=1.0):
z_hp = wo / z
p_hp = wo / p
if z_hp == float('inf'):
z_hp = []
# If lowpass had zeros at infinity, inverting moves them to origin.
z_hp = append(z_hp, zeros(degree))
z_hp = append(z_hp, np.zeros(degree))
# Cancel out gain change caused by inversion
k_hp = k * real(prod(-z) / prod(-p))
k_hp = k * (prod(-z) / prod(-p)).real
return z_hp, p_hp, k_hp
@ -1430,7 +1426,6 @@ def buttap(N):
k = 1
return z, p, k
def butter_bandpass(lowcut, highcut, fs, order=5):
nyq = 0.5 * fs
low = lowcut / nyq
@ -1438,7 +1433,6 @@ def butter_bandpass(lowcut, highcut, fs, order=5):
sos = butter(order, [low, high], analog=False, btype='band', output='sos')
return sos
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
sos = butter_bandpass(lowcut, highcut, fs, order=order)
y = spy.sosfilter(sos, data)
@ -1452,7 +1446,6 @@ def fft(x):
odd = fft(x[1::2])
return [even[m] + math.e**(-2j*math.pi*m/n)*odd[m] for m in range(n//2)] + [even[m] - math.e**(-2j*math.pi*m/n)*odd[m] for m in range(n//2)]
filter_dict = {'butter': [buttap, buttord],
'butterworth': [buttap, buttord]
}

View file

@ -4,3 +4,7 @@
5000
1
2.0
(array([0, 1, 2], dtype=int8),)
(array([1, 1, 1, 2, 2, 2], dtype=int8), array([0, 1, 2, 0, 1, 2], dtype=int8))
(array([1], dtype=int8),)
(array([3], dtype=int8),)

View file

@ -3,6 +3,8 @@ import sys
sys.path.append('.')
from snippets import numpy
from ulab import numpy as np
np.set_printoptions(threshold=100)
print (numpy.asarray([1]))
print (numpy.asarray([1.0, 2.0, 3j]))

View file

@ -1,3 +1,3 @@
array([1.0], dtype=float64)
array([1.0+0.0j, 2.0+0.0j, 0.0+3.0j], dtype=complex)
array([(4.0+0.0j, 3.0+0.0j, 1.0+0.0j, 2.0-2.0j, 2.0+2.0j, 2.0-1.0j, 2.0+1.0j, 2.0-1.0j, 2.0+1.0j, 1.0+1.0j, 1.0-1.0j], dtype=complex)
array([4.0+0.0j, 3.0+0.0j, 1.0+0.0j, 2.0-2.0j, 2.0+2.0j, 2.0-1.0j, 2.0+1.0j, 2.0-1.0j, 2.0+1.0j, 1.0+1.0j, 1.0-1.0j], dtype=complex)

View file

@ -1 +1,4 @@
array([[0, 0, 0],[0, 0, 0]])
array([[0.0, 0.0, 0.0],
[0.0, 0.0, 0.0]], dtype=float64)
array([[0.0+0.0j, 0.0+0.0j, 0.0+0.0j],
[0.0+0.0j, 0.0+0.0j, 0.0+0.0j]], dtype=complex)

View file

@ -0,0 +1,15 @@
import math
import sys
sys.path.append('.')
from ulab import numpy as np
from snippets import numpy
print(numpy.atleast_1d(1.0))
x = np.arange(9.0).reshape((3,3))
print(numpy.atleast_1d(x))
print(numpy.atleast_1d(x) is x)
print(numpy.atleast_1d(1, [3, 4]))

View file

@ -0,0 +1,6 @@
array([1.0], dtype=float64)
array([[0.0, 1.0, 2.0],
[3.0, 4.0, 5.0],
[6.0, 7.0, 8.0]], dtype=float64)
True
[array([1.0], dtype=float64), array([3.0, 4.0], dtype=float64)]

View file

@ -4,3 +4,8 @@ sys.path.append('.')
from snippets import numpy
from ulab import numpy as np
print(numpy.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]]))
print(numpy.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0))

View file

@ -0,0 +1,4 @@
array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0], dtype=float64)
array([[1.0, 2.0, 3.0],
[4.0, 5.0, 6.0],
[7.0, 8.0, 9.0]], dtype=float64)

View file

@ -0,0 +1,16 @@
import math
import sys
sys.path.append('.')
from snippets import numpy
from ulab import numpy as np
print(numpy.poly((0, 0, 0)))
print(numpy.poly((-1./2, 0, 1./2)))
print(numpy.poly((0.847, 0, 0.9883)))
#P = np.array([[0, 1./3], [-1./2, 0]])
#print(numpy.poly(P))

View file

@ -0,0 +1,3 @@
array([1.0, 0.0, 0.0, 0.0], dtype=float64)
array([1.0, 0.0, -0.25, 0.0], dtype=float64)
array([1.0, -1.8353, 0.8370901, 0.0], dtype=float64)

View file

@ -5,5 +5,13 @@ sys.path.append('.')
from snippets import numpy
from ulab import numpy as np
a = np.array([1, 2j, 3, 4j], dtype=np.complex)
print (numpy.isreal(a))
a = np.array([1+1j, 1+0j, 4.5, 3, 2, 2j], dtype=np.complex)
print(numpy.isreal(a))
a = np.array([1+2j, 2+1j], dtype=np.complex)
print(numpy.isreal(a))
print(numpy.isreal(1))
print(numpy.isreal(1j))

View file

@ -0,0 +1,4 @@
[False, True, True, True, True, False]
[False, False]
True
False

View file

@ -21,10 +21,39 @@ lowerCutoffFrequency_Hz = centerFrequency_Hz/math.sqrt(2)
upperCutoffFrequenc_Hz = centerFrequency_Hz*math.sqrt(2)
wn = np.array([ lowerCutoffFrequency_Hz, upperCutoffFrequenc_Hz])/nyquistRate
#print('butter: ', scipy.butter(N=4, Wn=wn, btype='bandpass', analog=False, output='ba'))
#print('butter: ', scipy.butter(N=4, Wn=wn, btype='bandpass', analog=False, output='zpk'))
print('butter: ', scipy.butter(N=4, Wn=wn, btype='bandpass', analog=False, output='sos'))
#print('butter: ', scipy.butter(N=4, Wn=wn, btype='bandstop', analog=False, output='ba'))
#print('butter: ', scipy.butter(N=4, Wn=wn, btype='lowpass', analog=False, output='ba'))
#print('butter: ', scipy.butter(N=4, Wn=wn, btype='highpass', analog=False, output='ba'))
z = []
p = np.array([-0.1564344650402309+0.9876883405951379j, -0.4539904997395468+0.8910065241883679j,
-0.7071067811865476+0.7071067811865475j, -0.8910065241883679+0.4539904997395467j, -0.9876883405951379+0.1564344650402309j,
-0.9876883405951379-0.1564344650402309j, -0.8910065241883679-0.4539904997395467j, -0.7071067811865476-0.7071067811865475j,
-0.4539904997395468-0.8910065241883679j, -0.1564344650402309-0.9876883405951379j], dtype=np.complex)
k = 1
wo = 0.1886352115099219
print(scipy.lp2hp_zpk(z,p,k,wo))
z = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], dtype=np.float)
p = np.array([-0.02950904840030544-0.1863127990340476j, -0.08563859394186457-0.1680752041469931j,
-0.1333852372292245-0.1333852372292244j, -0.1680752041469931-0.08563859394186453j, -0.1863127990340476-0.02950904840030543j,
-0.1863127990340476+0.02950904840030543j, -0.1680752041469931+0.08563859394186453j, -0.1333852372292245+0.1333852372292244j,
-0.08563859394186457+0.1680752041469931j, -0.02950904840030544+0.1863127990340476j], dtype=np.complex)
k = 1.0
fs = 2.0
print(scipy.bilinear_zpk(z,p,k,fs))
z = np.array([], dtype=np.float)
p = np.array([-0.3826834323650898+0.9238795325112868j,
-0.9238795325112868+0.3826834323650898j, -0.9238795325112868-0.3826834323650898j,
-0.3826834323650898-0.9238795325112868j], dtype=np.complex)
k = 1
wo = 0.03141673402115484
bw = 0.02221601345771878
print(scipy.lp2bs_zpk(z, p, k, wo=wo, bw=bw))
print(scipy.butter(N=4, Wn=wn, btype='bandpass', analog=False, output='ba'))
print(scipy.butter(N=4, Wn=wn, btype='bandpass', analog=False, output='zpk'))
print(scipy.butter(N=4, Wn=wn, btype='bandpass', analog=False, output='sos'))
print(scipy.butter(N=4, Wn=wn, btype='bandstop', analog=False, output='ba'))
print(scipy.butter(10, 15, 'lp', fs=1000, output='sos'))
print(scipy.butter(10, 15, 'hp', fs=1000, output='sos'))

View file

@ -1 +1,21 @@
butter: (array([9.375938694653579e-10, 0.0, -3.750375477861432e-09, 0.0, 5.625563216792147e-09, 0.0, -3.750375477861432e-09, 0.0, 9.375938694653579e-10], dtype=float64), array([1.0, -7.969992171296993, 27.79137077985833, -55.37836320535709, 68.97098091968704, -54.97798122059838, 27.39096409669159, -7.798371663621388, 0.9713924646368845], dtype=float64))
(array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], dtype=float64), array([-0.02950904840030542-0.1863127990340476j, -0.08563859394186453-0.1680752041469931j, -0.1333852372292245-0.1333852372292245j, -0.1680752041469931-0.08563859394186455j, -0.1863127990340476-0.02950904840030542j, -0.1863127990340476+0.02950904840030542j, -0.1680752041469931+0.08563859394186455j, -0.1333852372292245+0.1333852372292245j, -0.08563859394186453+0.1680752041469931j, -0.02950904840030542+0.1863127990340476j], dtype=complex), 1.0)
(array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=float64), array([0.9811181553845975-0.09160115148355542j, 0.9547700993580936-0.08041542979283999j, 0.9334461398558574-0.06239272587315813j, 0.918541250895572-0.03941895649644052j, 0.9108945995807357-0.01346977255018339j, 0.9108945995807357+0.01346977255018339j, 0.918541250895572+0.03941895649644052j, 0.9334461398558574+0.06239272587315813j, 0.9547700993580936+0.08041542979283999j, 0.9811181553845975+0.09160115148355542j], dtype=complex), 0.73979400584056)
(array([0.0+0.03141673402115484j, 0.0+0.03141673402115484j, 0.0+0.03141673402115484j, 0.0+0.03141673402115484j, 0.0-0.03141673402115484j, 0.0-0.03141673402115484j, 0.0-0.03141673402115484j, 0.0-0.03141673402115484j], dtype=complex), array([-0.002920960939069395+0.02254040782120927j, -0.008809831448862637+0.02578034941498454j, -0.008809831448862637-0.02578034941498454j, -0.002920960939069395-0.02254040782120927j, -0.005580739344399452-0.04306532794879095j, -0.01171508867871904-0.03428204969845339j, -0.01171508867871904+0.03428204969845339j, -0.005580739344399452+0.04306532794879095j], dtype=complex), 1.0)
(array([9.375938694653579e-10, 0.0, -3.750375477861432e-09, 0.0, 5.625563216792147e-09, 0.0, -3.750375477861432e-09, 0.0, 9.375938694653579e-10], dtype=float64), array([1.0, -7.969992171296993, 27.79137077985833, -55.37836320535709, 68.97098091968704, -54.97798122059838, 27.39096409669159, -7.798371663621388, 0.9713924646368845], dtype=float64))
(array([1.0+0.0j, 1.0+0.0j, 1.0+0.0j, 1.0+0.0j, -1.0+0.0j, -1.0+0.0j, -1.0+0.0j, -1.0+0.0j], dtype=complex), array([0.9984772174419884-0.01125340518638924j, 0.9955222362276896-0.01283305087503436j, 0.9955222362276896+0.01283305087503436j, 0.9984772174419884+0.01125340518638924j, 0.9969826844852829+0.02147022362342683j, 0.9940139474935357+0.01703981557421542j, 0.9940139474935357-0.01703981557421542j, 0.9969826844852829-0.02147022362342683j], dtype=complex), 9.375938694653579e-10)
array([[9.375938694653579e-10, 1.875187738930716e-09, 9.375938694653579e-10, 1.0, -1.988027894987071, 0.9883540831264847],
[1.0, 2.0, 1.0, 1.0, -1.991044472455379, 0.9912292100185409],
[1.0, -2.0, 1.0, 1.0, -1.993965368970566, 0.9944354436659212],
[1.0, -2.0, 1.0, 1.0, -1.996954434883977, 0.9970833928789848]], dtype=float64)
(array([0.9855924434759455, -7.883766817056333, 27.59075239283292, -55.17858731338059, 68.97201858825612, -55.17858731338059, 27.59075239283291, -7.883766817056331, 0.9855924434759455], dtype=float64), array([1.0, -7.969992171296993, 27.79137077985833, -55.37836320535709, 68.97098091968704, -54.97798122059838, 27.39096409669159, -7.798371663621388, 0.9713924646368845], dtype=float64))
array([[1.609898589131747e-13, 3.219797178263494e-13, 1.609898589131747e-13, 1.0, 0.0, 0.0],
[1.0, 2.0, 1.0, 1.0, -1.821789199161471, 0.8299104063179026],
[1.0, 2.0, 1.0, 1.0, -1.837082501791144, 0.8452718837280704],
[1.0, 2.0, 1.0, 1.0, -1.866892279711715, 0.8752145482536838],
[1.0, 2.0, 1.0, 1.0, -1.909540198716187, 0.918052583977031],
[1.0, -1.0, 0.0, 1.0, -1.962236310769195, 0.9709836057783884]], dtype=float64)
array([[0.73979400584056, -1.47958801168112, 0.73979400584056, 1.0, -1.821789199161471, 0.8299104063179026],
[1.0, -2.0, 1.0, 1.0, -1.837082501791144, 0.8452718837280704],
[1.0, -2.0, 1.0, 1.0, -1.866892279711715, 0.8752145482536838],
[1.0, -2.0, 1.0, 1.0, -1.909540198716187, 0.918052583977031],
[1.0, -2.0, 1.0, 1.0, -1.962236310769195, 0.9709836057783884]], dtype=float64)

0
test-common.sh Normal file → Executable file
View file

18
test-snippets.sh Executable file
View file

@ -0,0 +1,18 @@
#!/bin/sh
set -e
micropython="$1"
for level1 in numpy scipy;
do
for level2 in core lib signal; do
rm -f *.exp
if ! env MICROPY_MICROPYTHON="$micropython" ./run-tests -d snippets/tests/"$level1"/"$level2"; then
for exp in *.exp; do
testbase=$(basename $exp .exp);
echo -e "\nFAILURE $testbase";
diff -u $testbase.exp $testbase.out;
done
exit 1
fi
done
done

View file

@ -1,32 +0,0 @@
import math
from ulab import numpy
from ulab import numpy as np
from ulab import scipy as spy
def to_tuple(a):
result = []
for x in a:
result.append([x.real, x.imag])
return result
def to_complex(a):
result = np.array([], dtype=np.complex)
for x in a:
t = np.array([complex(x[0] + x[1] * 1j)], dtype=np.complex)
result = np.concatenate((result, t), axis=0)
return result
roots = np.array([0.99847722-0.01125341j, 0.99552224-0.01283305j, 0.99552224+0.01283305j, 0.99847722+0.01125341j,
0.99698268+0.02147022j, 0.9940139500000001+0.01703982j, 0.9940139500000001-0.01703982j, 0.99698268-0.02147022j], dtype=np.complex)
print(roots)
r = to_tuple(roots)
s = sorted(r,key=lambda x:(x[0], x[1]))
f = to_complex(s)
print(f)
#print()
#print(sorted(r,key=lambda x:(x[0], x[1])))
#print(np.sort_complex(np.array([-0.9800000000000001+1.0j, 0.9800000000000001-1.0j, 0.99+1.0j, 0.99-1.0j], dtype=np.complex)))

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

View file

@ -1,70 +0,0 @@
import math
from ulab import numpy
from ulab import numpy as np
from ulab import scipy as spy
def asarray(a, dtype=None):
if isinstance(a,(np.ndarray)):
return a
return a
def atleast_1d(*arys):
res = []
for ary in arys:
ary = asarray(ary)
if not isinstance(ary,(np.ndarray)):
ary = np.array([ary])
result = ary.reshape((1,))
else:
result = ary
res.append(result)
if len(res) == 1:
return res[0]
else:
return res
def poly(seq_of_zeros):
seq_of_zeros = atleast_1d(seq_of_zeros)
sh = seq_of_zeros.shape
if len(sh) == 2 and sh[0] == sh[1] and sh[0] != 0:
seq_of_zeros = eigvals(seq_of_zeros)
elif len(sh) == 1:
dt = seq_of_zeros.dtype
# Let object arrays slip through, e.g. for arbitrary precision
if dt != object:
seq_of_zeros = seq_of_zeros #seq_of_zeros.astype(mintypecode(dt.char))
else:
raise ValueError("input must be 1d or non-empty square 2d array.")
if len(seq_of_zeros) == 0:
return 1.0
dt = seq_of_zeros.dtype
print(dt)
a = np.ones((1,), dtype=dt)
print(a)
print(seq_of_zeros)
for k in range(len(seq_of_zeros)):
a = np.convolve(a, np.array([1, -seq_of_zeros[k]], dtype=dt))
print(a)
#if issubclass(a.dtype.type, NX.complexfloating):
# if complex roots are all complex conjugates, the roots are real.
roots = asarray(seq_of_zeros, complex)
print('roots',roots)
p = np.sort_complex(roots)
print(p)
q = np.conjugate(roots)
s = np.sort_complex(q)
print(s)
#if NX.all(NX.sort(roots) == NX.sort(roots.conjugate())):
a = a.real.copy()
return a
p = np.array([-0.3826834323650898+0.9238795325112868j, -0.9238795325112868+0.3826834323650898j,
-0.9238795325112868-0.3826834323650898j, -0.3826834323650898-0.9238795325112868j], dtype=np.complex)
print(poly(p))