From e96fa5b305ffb0acc92e00eff676ad74325f2a5c Mon Sep 17 00:00:00 2001 From: Warren Weckesser Date: Sun, 29 May 2011 23:42:28 -0500 Subject: [PATCH] STY: signal: big PEP8 cleanup, plus a BUG FIX: there was a comma missing in the list of coefficients for N=11 in the besselap function in filter_design.py. --- scipy/signal/__init__.py | 2 +- scipy/signal/filter_design.py | 1262 +++++++++++++++-------------- scipy/signal/fir_filter_design.py | 54 +- scipy/signal/info.py | 14 +- scipy/signal/ltisys.py | 111 +-- scipy/signal/setup.py | 14 +- scipy/signal/setupscons.py | 4 +- scipy/signal/signaltools.py | 4 +- scipy/signal/waveforms.py | 11 +- scipy/signal/wavelets.py | 70 +- scipy/signal/windows.py | 85 +- 11 files changed, 864 insertions(+), 767 deletions(-) diff --git a/scipy/signal/__init__.py b/scipy/signal/__init__.py index a0eaf790242a..13e810afa620 100644 --- a/scipy/signal/__init__.py +++ b/scipy/signal/__init__.py @@ -20,6 +20,6 @@ from spectral import * from wavelets import * -__all__ = filter(lambda s:not s.startswith('_'),dir()) +__all__ = filter(lambda s: not s.startswith('_'), dir()) from numpy.testing import Tester test = Tester().test diff --git a/scipy/signal/filter_design.py b/scipy/signal/filter_design.py index f69a7ed9a6da..42cb0aca3b1f 100644 --- a/scipy/signal/filter_design.py +++ b/scipy/signal/filter_design.py @@ -25,22 +25,28 @@ class BadCoefficients(UserWarning): abs = absolute + def findfreqs(num, den, N): - ep = atleast_1d(roots(den))+0j - tz = atleast_1d(roots(num))+0j + ep = atleast_1d(roots(den)) + 0j + tz = atleast_1d(roots(num)) + 0j if len(ep) == 0: - ep = atleast_1d(-1000)+0j + ep = atleast_1d(-1000) + 0j - ez = r_['-1',numpy.compress(ep.imag >=0, ep,axis=-1), numpy.compress((abs(tz) < 1e5) & (tz.imag >=0),tz,axis=-1)] + ez = r_['-1', + numpy.compress(ep.imag >= 0, ep, axis=-1), + numpy.compress((abs(tz) < 1e5) & (tz.imag >= 0), tz, axis=-1)] integ = abs(ez) < 1e-10 - hfreq = numpy.around(numpy.log10(numpy.max(3*abs(ez.real + integ)+1.5*ez.imag))+0.5) - lfreq = numpy.around(numpy.log10(0.1*numpy.min(abs(real(ez+integ))+2*ez.imag))-0.5) + hfreq = numpy.around(numpy.log10(numpy.max(3 * abs(ez.real + integ) + + 1.5 * ez.imag)) + 0.5) + lfreq = numpy.around(numpy.log10(0.1 * numpy.min(abs(real(ez + integ)) + + 2 * ez.imag)) - 0.5) w = logspace(lfreq, hfreq, N) return w + def freqs(b, a, worN=None, plot=None): """ Compute frequency response of analog filter. @@ -87,19 +93,20 @@ def freqs(b, a, worN=None, plot=None): """ if worN is None: - w = findfreqs(b,a,200) + w = findfreqs(b, a, 200) elif isinstance(worN, types.IntType): N = worN - w = findfreqs(b,a,N) + w = findfreqs(b, a, N) else: w = worN w = atleast_1d(w) - s = 1j*w + s = 1j * w h = polyval(b, s) / polyval(a, s) if not plot is None: plot(w, h) return w, h + def freqz(b, a=1, worN=None, whole=0, plot=None): """ Compute the frequency response of a digital filter. @@ -168,9 +175,9 @@ def freqz(b, a=1, worN=None, whole=0, plot=None): >>> plt.show() """ - b, a = map(atleast_1d, (b,a)) + b, a = map(atleast_1d, (b, a)) if whole: - lastpoint = 2*pi + lastpoint = 2 * pi else: lastpoint = pi if worN is None: @@ -182,12 +189,13 @@ def freqz(b, a=1, worN=None, whole=0, plot=None): else: w = worN w = atleast_1d(w) - zm1 = exp(-1j*w) + zm1 = exp(-1j * w) h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1) if not plot is None: plot(w, h) return w, h + def tf2zpk(b, a): """Return zero, pole, gain (z,p,k) representation from a numerator, denominator representation of a linear filter. @@ -211,15 +219,16 @@ def tf2zpk(b, a): If some values of b are too close to 0, they are removed. In that case, a BadCoefficients warning is emitted. """ - b,a = normalize(b,a) - b = (b+0.0) / a[0] - a = (a+0.0) / a[0] + b, a = normalize(b, a) + b = (b + 0.0) / a[0] + a = (a + 0.0) / a[0] k = b[0] b /= b[0] z = roots(b) p = roots(a) return z, p, k + def zpk2tf(z, p, k): """Return polynomial transfer function representation from zeros and poles @@ -245,9 +254,9 @@ def zpk2tf(z, p, k): k = atleast_1d(k) if len(z.shape) > 1: temp = poly(z[0]) - b = zeros((z.shape[0], z.shape[1]+1), temp.dtype.char) + b = zeros((z.shape[0], z.shape[1] + 1), temp.dtype.char) if len(k) == 1: - k = [k[0]]*z.shape[0] + k = [k[0]] * z.shape[0] for i in range(z.shape[0]): b[i] = k[i] * poly(z[i]) else: @@ -255,28 +264,30 @@ def zpk2tf(z, p, k): a = atleast_1d(poly(p)) return b, a + def normalize(b, a): """Normalize polynomial representation of a transfer function. If values of b are too close to 0, they are removed. In that case, a BadCoefficients warning is emitted. """ - b,a = map(atleast_1d,(b,a)) + b, a = map(atleast_1d, (b, a)) if len(a.shape) != 1: raise ValueError("Denominator polynomial must be rank-1 array.") if len(b.shape) > 2: - raise ValueError("Numerator polynomial must be rank-1 or rank-2 array.") + raise ValueError("Numerator polynomial must be rank-1 or" + " rank-2 array.") if len(b.shape) == 1: - b = asarray([b],b.dtype.char) + b = asarray([b], b.dtype.char) while a[0] == 0.0 and len(a) > 1: a = a[1:] outb = b * (1.0) / a[0] outa = a * (1.0) / a[0] - if allclose(outb[:,0], 0, rtol=1e-14): + if allclose(outb[:, 0], 0, rtol=1e-14): warnings.warn("Badly conditioned filter coefficients (numerator): the " "results may be meaningless", BadCoefficients) - while allclose(outb[:,0], 0, rtol=1e-14) and (outb.shape[-1] > 1): - outb = outb[:,1:] + while allclose(outb[:, 0], 0, rtol=1e-14) and (outb.shape[-1] > 1): + outb = outb[:, 1:] if outb.shape[0] == 1: outb = outb[0] return outb, outa @@ -286,26 +297,27 @@ def lp2lp(b, a, wo=1.0): """Return a low-pass filter with cutoff frequency `wo` from a low-pass filter prototype with unity cutoff frequency. """ - a,b = map(atleast_1d,(a,b)) + a, b = map(atleast_1d, (a, b)) try: wo = float(wo) except TypeError: wo = float(wo[0]) d = len(a) n = len(b) - M = max((d,n)) - pwo = pow(wo,numpy.arange(M-1,-1,-1)) - start1 = max((n-d,0)) - start2 = max((d-n,0)) - b = b * pwo[start1]/pwo[start2:] - a = a * pwo[start1]/pwo[start1:] + M = max((d, n)) + pwo = pow(wo, numpy.arange(M - 1, -1, -1)) + start1 = max((n - d, 0)) + start2 = max((d - n, 0)) + b = b * pwo[start1] / pwo[start2:] + a = a * pwo[start1] / pwo[start1:] return normalize(b, a) + def lp2hp(b, a, wo=1.0): """Return a high-pass filter with cutoff frequency `wo` from a low-pass filter prototype with unity cutoff frequency. """ - a,b = map(atleast_1d,(a,b)) + a, b = map(atleast_1d, (a, b)) try: wo = float(wo) except TypeError: @@ -313,118 +325,126 @@ def lp2hp(b, a, wo=1.0): d = len(a) n = len(b) if wo != 1: - pwo = pow(wo,numpy.arange(max((d,n)))) + pwo = pow(wo, numpy.arange(max((d, n)))) else: - pwo = numpy.ones(max((d,n)),b.dtype.char) + pwo = numpy.ones(max((d, n)), b.dtype.char) if d >= n: outa = a[::-1] * pwo - outb = resize(b,(d,)) + outb = resize(b, (d,)) outb[n:] = 0.0 outb[:n] = b[::-1] * pwo[:n] else: outb = b[::-1] * pwo - outa = resize(a,(n,)) + outa = resize(a, (n,)) outa[d:] = 0.0 outa[:d] = a[::-1] * pwo[:d] return normalize(outb, outa) + def lp2bp(b, a, wo=1.0, bw=1.0): """Return a band-pass filter with center frequency `wo` and bandwidth `bw` from a low-pass filter prototype with unity cutoff frequency. """ - a,b = map(atleast_1d,(a,b)) + a, b = map(atleast_1d, (a, b)) D = len(a) - 1 N = len(b) - 1 - artype = mintypecode((a,b)) - ma = max([N,D]) + artype = mintypecode((a, b)) + ma = max([N, D]) Np = N + ma Dp = D + ma - bprime = numpy.zeros(Np+1,artype) - aprime = numpy.zeros(Dp+1,artype) - wosq = wo*wo - for j in range(Np+1): + bprime = numpy.zeros(Np + 1, artype) + aprime = numpy.zeros(Dp + 1, artype) + wosq = wo * wo + for j in range(Np + 1): val = 0.0 - for i in range(0,N+1): - for k in range(0,i+1): - if ma-i+2*k == j: - val += comb(i,k)*b[N-i]*(wosq)**(i-k) / bw**i - bprime[Np-j] = val - for j in range(Dp+1): + for i in range(0, N + 1): + for k in range(0, i + 1): + if ma - i + 2 * k == j: + val += comb(i, k) * b[N - i] * (wosq) ** (i - k) / bw ** i + bprime[Np - j] = val + for j in range(Dp + 1): val = 0.0 - for i in range(0,D+1): - for k in range(0,i+1): - if ma-i+2*k == j: - val += comb(i,k)*a[D-i]*(wosq)**(i-k) / bw**i - aprime[Dp-j] = val + for i in range(0, D + 1): + for k in range(0, i + 1): + if ma - i + 2 * k == j: + val += comb(i, k) * a[D - i] * (wosq) ** (i - k) / bw ** i + aprime[Dp - j] = val return normalize(bprime, aprime) + def lp2bs(b, a, wo=1, bw=1): """Return a band-stop filter with center frequency `wo` and bandwidth `bw` from a low-pass filter prototype with unity cutoff frequency. """ - a,b = map(atleast_1d,(a,b)) + a, b = map(atleast_1d, (a, b)) D = len(a) - 1 N = len(b) - 1 - artype = mintypecode((a,b)) - M = max([N,D]) + artype = mintypecode((a, b)) + M = max([N, D]) Np = M + M Dp = M + M - bprime = numpy.zeros(Np+1,artype) - aprime = numpy.zeros(Dp+1,artype) - wosq = wo*wo - for j in range(Np+1): + bprime = numpy.zeros(Np + 1, artype) + aprime = numpy.zeros(Dp + 1, artype) + wosq = wo * wo + for j in range(Np + 1): val = 0.0 - for i in range(0,N+1): - for k in range(0,M-i+1): - if i+2*k == j: - val += comb(M-i,k)*b[N-i]*(wosq)**(M-i-k) * bw**i - bprime[Np-j] = val - for j in range(Dp+1): + for i in range(0, N + 1): + for k in range(0, M - i + 1): + if i + 2 * k == j: + val += (comb(M - i, k) * b[N - i] * + (wosq) ** (M - i - k) * bw ** i) + bprime[Np - j] = val + for j in range(Dp + 1): val = 0.0 - for i in range(0,D+1): - for k in range(0,M-i+1): - if i+2*k == j: - val += comb(M-i,k)*a[D-i]*(wosq)**(M-i-k) * bw**i - aprime[Dp-j] = val + for i in range(0, D + 1): + for k in range(0, M - i + 1): + if i + 2 * k == j: + val += (comb(M - i, k) * a[D - i] * + (wosq) ** (M - i - k) * bw ** i) + aprime[Dp - j] = val return normalize(bprime, aprime) + def bilinear(b, a, fs=1.0): - """Return a digital filter from an analog filter using the bilinear transform. + """Return a digital filter from an analog one using a bilinear transform. The bilinear transform substitutes ``(z-1) / (z+1``) for ``s``. """ - fs =float(fs) - a,b = map(atleast_1d,(a,b)) + fs = float(fs) + a, b = map(atleast_1d, (a, b)) D = len(a) - 1 N = len(b) - 1 artype = float - M = max([N,D]) + M = max([N, D]) Np = M Dp = M - bprime = numpy.zeros(Np+1,artype) - aprime = numpy.zeros(Dp+1,artype) - for j in range(Np+1): + bprime = numpy.zeros(Np + 1, artype) + aprime = numpy.zeros(Dp + 1, artype) + for j in range(Np + 1): val = 0.0 - for i in range(N+1): - for k in range(i+1): - for l in range(M-i+1): - if k+l == j: - val += comb(i,k)*comb(M-i,l)*b[N-i]*pow(2*fs,i)*(-1)**k + for i in range(N + 1): + for k in range(i + 1): + for l in range(M - i + 1): + if k + l == j: + val += (comb(i, k) * comb(M - i, l) * b[N - i] * + pow(2 * fs, i) * (-1) ** k) bprime[j] = real(val) - for j in range(Dp+1): + for j in range(Dp + 1): val = 0.0 - for i in range(D+1): - for k in range(i+1): - for l in range(M-i+1): - if k+l == j: - val += comb(i,k)*comb(M-i,l)*a[D-i]*pow(2*fs,i)*(-1)**k + for i in range(D + 1): + for k in range(i + 1): + for l in range(M - i + 1): + if k + l == j: + val += (comb(i, k) * comb(M - i, l) * a[D - i] * + pow(2 * fs, i) * (-1) ** k) aprime[j] = real(val) return normalize(bprime, aprime) + def iirdesign(wp, ws, gpass, gstop, analog=0, ftype='ellip', output='ba'): """Complete IIR digital and analog filter design. @@ -478,22 +498,26 @@ def iirdesign(wp, ws, gpass, gstop, analog=0, ftype='ellip', output='ba'): except KeyError: raise ValueError("Invalid IIR filter type: %s" % ftype) except IndexError: - raise ValueError("%s does not have order selection use iirfilter function." % ftype) + raise ValueError(("%s does not have order selection use " + "iirfilter function.") % ftype) wp = atleast_1d(wp) ws = atleast_1d(ws) - band_type = 2*(len(wp)-1) - band_type +=1 + band_type = 2 * (len(wp) - 1) + band_type += 1 if wp[0] >= ws[0]: band_type += 1 - btype = {1:'lowpass', 2:'highpass', 3:'bandstop', 4:'bandpass'}[band_type] + btype = {1: 'lowpass', 2: 'highpass', + 3: 'bandstop', 4: 'bandpass'}[band_type] N, Wn = ordfunc(wp, ws, gpass, gstop, analog=analog) - return iirfilter(N, Wn, rp=gpass, rs=gstop, analog=analog, btype=btype, ftype=ftype, output=output) + return iirfilter(N, Wn, rp=gpass, rs=gstop, analog=analog, btype=btype, + ftype=ftype, output=output) -def iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=0, ftype='butter', output='ba'): +def iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=0, + ftype='butter', output='ba'): """IIR digital and analog filter design given order and critical points. Design an Nth order lowpass digital or analog filter and return the filter @@ -550,10 +574,10 @@ def iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=0, ftype='butter', o if output not in ['ba', 'zpk']: raise ValueError("%s is not a valid output form." % output) - #pre-warp frequencies for digital filter design + # pre-warp frequencies for digital filter design if not analog: fs = 2.0 - warped = 2*fs*tan(pi*Wn/fs) + warped = 2 * fs * tan(pi * Wn / fs) else: warped = Wn @@ -562,36 +586,38 @@ def iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=0, ftype='butter', o wo = warped else: bw = warped[1] - warped[0] - wo = sqrt(warped[0]*warped[1]) + wo = sqrt(warped[0] * warped[1]) # Get analog lowpass prototype if typefunc in [buttap, besselap]: z, p, k = typefunc(N) elif typefunc == cheb1ap: if rp is None: - raise ValueError("passband ripple (rp) must be provided to design a Chebyshev I filter.") + raise ValueError("passband ripple (rp) must be provided to " + "design a Chebyshev I filter.") z, p, k = typefunc(N, rp) elif typefunc == cheb2ap: if rs is None: - raise ValueError("stopband atteunatuion (rs) must be provided to design an Chebyshev II filter.") + raise ValueError("stopband atteunatuion (rs) must be provided to " + "design an Chebyshev II filter.") z, p, k = typefunc(N, rs) else: # Elliptic filters if rs is None or rp is None: - raise ValueError("Both rp and rs must be provided to design an elliptic filter.") + raise ValueError("Both rp and rs must be provided to design an " + "elliptic filter.") z, p, k = typefunc(N, rp, rs) - b, a = zpk2tf(z,p,k) + b, a = zpk2tf(z, p, k) # transform to lowpass, bandpass, highpass, or bandstop if btype == 'lowpass': - b, a = lp2lp(b,a,wo=wo) + b, a = lp2lp(b, a, wo=wo) elif btype == 'highpass': - b, a = lp2hp(b,a,wo=wo) + b, a = lp2hp(b, a, wo=wo) elif btype == 'bandpass': - b, a = lp2bp(b,a,wo=wo,bw=bw) - else: # 'bandstop' - b, a = lp2bs(b,a,wo=wo,bw=bw) - + b, a = lp2bp(b, a, wo=wo, bw=bw) + else: # 'bandstop' + b, a = lp2bs(b, a, wo=wo, bw=bw) # Find discrete equivalent if necessary if not analog: @@ -599,9 +625,9 @@ def iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=0, ftype='butter', o # Transform to proper out type (pole-zero, state-space, numer-denom) if output == 'zpk': - return tf2zpk(b,a) + return tf2zpk(b, a) else: - return b,a + return b, a def butter(N, Wn, btype='low', analog=0, output='ba'): @@ -614,7 +640,9 @@ def butter(N, Wn, btype='low', analog=0, output='ba'): -------- buttord. """ - return iirfilter(N, Wn, btype=btype, analog=analog, output=output, ftype='butter') + return iirfilter(N, Wn, btype=btype, analog=analog, + output=output, ftype='butter') + def cheby1(N, rp, Wn, btype='low', analog=0, output='ba'): """Chebyshev type I digital and analog filter design. @@ -626,7 +654,9 @@ def cheby1(N, rp, Wn, btype='low', analog=0, output='ba'): -------- cheb1ord. """ - return iirfilter(N, Wn, rp=rp, btype=btype, analog=analog, output=output, ftype='cheby1') + return iirfilter(N, Wn, rp=rp, btype=btype, analog=analog, + output=output, ftype='cheby1') + def cheby2(N, rs, Wn, btype='low', analog=0, output='ba'): """Chebyshev type I digital and analog filter design. @@ -638,7 +668,9 @@ def cheby2(N, rs, Wn, btype='low', analog=0, output='ba'): -------- cheb2ord. """ - return iirfilter(N, Wn, rs=rs, btype=btype, analog=analog, output=output, ftype='cheby2') + return iirfilter(N, Wn, rs=rs, btype=btype, analog=analog, + output=output, ftype='cheby2') + def ellip(N, rp, rs, Wn, btype='low', analog=0, output='ba'): """Elliptic (Cauer) digital and analog filter design. @@ -650,7 +682,9 @@ def ellip(N, rp, rs, Wn, btype='low', analog=0, output='ba'): -------- ellipord. """ - return iirfilter(N, Wn, rs=rs, rp=rp, btype=btype, analog=analog, output=output, ftype='elliptic') + return iirfilter(N, Wn, rs=rs, rp=rp, btype=btype, analog=analog, + output=output, ftype='elliptic') + def bessel(N, Wn, btype='low', analog=0, output='ba'): """Bessel digital and analog filter design. @@ -659,12 +693,14 @@ def bessel(N, Wn, btype='low', analog=0, output='ba'): filter coefficients in (B,A) or (Z,P,K) form. """ - return iirfilter(N, Wn, btype=btype, analog=analog, output=output, ftype='bessel') + return iirfilter(N, Wn, btype=btype, analog=analog, + output=output, ftype='bessel') def maxflat(): pass + def yulewalk(): pass @@ -699,29 +735,31 @@ def band_stop_obj(wp, ind, passb, stopb, gpass, gstop, type): passbC = passb.copy() passbC[ind] = wp - nat = stopb*(passbC[0]-passbC[1]) / (stopb**2 - passbC[0]*passbC[1]) + nat = (stopb * (passbC[0] - passbC[1]) / + (stopb ** 2 - passbC[0] * passbC[1])) nat = min(abs(nat)) if type == 'butter': - GSTOP = 10**(0.1*abs(gstop)) - GPASS = 10**(0.1*abs(gpass)) - n = (log10((GSTOP-1.0)/(GPASS-1.0)) / (2*log10(nat))) + GSTOP = 10 ** (0.1 * abs(gstop)) + GPASS = 10 ** (0.1 * abs(gpass)) + n = (log10((GSTOP - 1.0) / (GPASS - 1.0)) / (2 * log10(nat))) elif type == 'cheby': - GSTOP = 10**(0.1*abs(gstop)) - GPASS = 10**(0.1*abs(gpass)) - n = arccosh(sqrt((GSTOP-1.0)/(GPASS-1.0))) / arccosh(nat) + GSTOP = 10 ** (0.1 * abs(gstop)) + GPASS = 10 ** (0.1 * abs(gpass)) + n = arccosh(sqrt((GSTOP - 1.0) / (GPASS - 1.0))) / arccosh(nat) elif type == 'ellip': - GSTOP = 10**(0.1*gstop) - GPASS = 10**(0.1*gpass) - arg1 = sqrt( (GPASS-1.0) / (GSTOP-1.0) ) + GSTOP = 10 ** (0.1 * gstop) + GPASS = 10 ** (0.1 * gpass) + arg1 = sqrt((GPASS - 1.0) / (GSTOP - 1.0)) arg0 = 1.0 / nat - d0 = special.ellipk([arg0**2, 1-arg0**2]) - d1 = special.ellipk([arg1**2, 1-arg1**2]) - n = (d0[0]*d1[1] / (d0[1]*d1[0])) + d0 = special.ellipk([arg0 ** 2, 1 - arg0 ** 2]) + d1 = special.ellipk([arg1 ** 2, 1 - arg1 ** 2]) + n = (d0[0] * d1[1] / (d0[1] * d1[0])) else: raise ValueError("Incorrect type: %s" % type) return n + def buttord(wp, ws, gpass, gstop, analog=0): """Butterworth filter order selection. @@ -760,46 +798,50 @@ def buttord(wp, ws, gpass, gstop, analog=0): wp = atleast_1d(wp) ws = atleast_1d(ws) - filter_type = 2*(len(wp)-1) - filter_type +=1 + filter_type = 2 * (len(wp) - 1) + filter_type += 1 if wp[0] >= ws[0]: filter_type += 1 # Pre-warp frequencies if not analog: - passb = tan(wp*pi/2.0) - stopb = tan(ws*pi/2.0) + passb = tan(wp * pi / 2.0) + stopb = tan(ws * pi / 2.0) else: - passb = wp*1.0 - stopb = ws*1.0 + passb = wp * 1.0 + stopb = ws * 1.0 if filter_type == 1: # low nat = stopb / passb elif filter_type == 2: # high nat = passb / stopb elif filter_type == 3: # stop - wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0]-1e-12, - args=(0,passb,stopb,gpass,gstop,'butter'), + wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0] - 1e-12, + args=(0, passb, stopb, gpass, gstop, + 'butter'), disp=0) passb[0] = wp0 - wp1 = optimize.fminbound(band_stop_obj, stopb[1]+1e-12, passb[1], - args=(1,passb,stopb,gpass,gstop,'butter'), + wp1 = optimize.fminbound(band_stop_obj, stopb[1] + 1e-12, passb[1], + args=(1, passb, stopb, gpass, gstop, + 'butter'), disp=0) passb[1] = wp1 - nat = (stopb * (passb[0]-passb[1])) / (stopb**2 - passb[0]*passb[1]) + nat = ((stopb * (passb[0] - passb[1])) / + (stopb ** 2 - passb[0] * passb[1])) elif filter_type == 4: # pass - nat = (stopb**2 - passb[0]*passb[1]) / (stopb* (passb[0]-passb[1])) + nat = ((stopb ** 2 - passb[0] * passb[1]) / + (stopb * (passb[0] - passb[1]))) nat = min(abs(nat)) - GSTOP = 10**(0.1*abs(gstop)) - GPASS = 10**(0.1*abs(gpass)) - ord = int(ceil( log10((GSTOP-1.0)/(GPASS-1.0)) / (2*log10(nat)))) + GSTOP = 10 ** (0.1 * abs(gstop)) + GPASS = 10 ** (0.1 * abs(gpass)) + ord = int(ceil(log10((GSTOP - 1.0) / (GPASS - 1.0)) / (2 * log10(nat)))) # Find the butterworth natural frequency W0 (or the "3dB" frequency") # to give exactly gstop at nat. W0 will be between 1 and nat try: - W0 = nat / ( ( 10**(0.1*abs(gstop))-1)**(1.0/(2.0*ord))) + W0 = nat / ((10 ** (0.1 * abs(gstop)) - 1) ** (1.0 / (2.0 * ord))) except ZeroDivisionError: W0 = nat print "Warning, order is zero...check input parametegstop." @@ -808,27 +850,27 @@ def buttord(wp, ws, gpass, gstop, analog=0): # to the original analog filter if filter_type == 1: # low - WN = W0*passb - elif filter_type == 2: # high + WN = W0 * passb + elif filter_type == 2: # high WN = passb / W0 elif filter_type == 3: # stop - WN = numpy.zeros(2,float) - WN[0] = ((passb[1] - passb[0]) + sqrt((passb[1] - passb[0])**2 + \ - 4*W0**2 * passb[0] * passb[1])) / (2*W0) - WN[1] = ((passb[1] - passb[0]) - sqrt((passb[1] - passb[0])**2 + \ - 4*W0**2 * passb[0] * passb[1])) / (2*W0) + WN = numpy.zeros(2, float) + discr = sqrt((passb[1] - passb[0]) ** 2 + + 4 * W0 ** 2 * passb[0] * passb[1]) + WN[0] = ((passb[1] - passb[0]) + discr) / (2 * W0) + WN[1] = ((passb[1] - passb[0]) - discr) / (2 * W0) WN = numpy.sort(abs(WN)) - elif filter_type == 4: # pass - W0 = numpy.array([-W0, W0],float) - WN = -W0 * (passb[1]-passb[0]) / 2.0 + sqrt(W0**2 / 4.0 * \ - (passb[1]-passb[0])**2 + \ - passb[0]*passb[1]) + elif filter_type == 4: # pass + W0 = numpy.array([-W0, W0], float) + WN = (-W0 * (passb[1] - passb[0]) / 2.0 + + sqrt(W0 ** 2 / 4.0 * (passb[1] - passb[0]) ** 2 + + passb[0] * passb[1])) WN = numpy.sort(abs(WN)) else: raise ValueError("Bad type: %s" % filter_type) if not analog: - wn = (2.0/pi)*arctan(WN) + wn = (2.0 / pi) * arctan(WN) else: wn = WN @@ -836,6 +878,7 @@ def buttord(wp, ws, gpass, gstop, analog=0): wn = wn[0] return ord, wn + def cheb1ord(wp, ws, gpass, gstop, analog=0): """Chebyshev type I filter order selection. @@ -873,7 +916,7 @@ def cheb1ord(wp, ws, gpass, gstop, analog=0): """ wp = atleast_1d(wp) ws = atleast_1d(ws) - filter_type = 2*(len(wp)-1) + filter_type = 2 * (len(wp) - 1) if wp[0] < ws[0]: filter_type += 1 else: @@ -881,36 +924,41 @@ def cheb1ord(wp, ws, gpass, gstop, analog=0): # Pre-wagpass frequencies if not analog: - passb = tan(pi*wp/2.) - stopb = tan(pi*ws/2.) + passb = tan(pi * wp / 2.) + stopb = tan(pi * ws / 2.) else: - passb = wp*1.0 - stopb = ws*1.0 + passb = wp * 1.0 + stopb = ws * 1.0 if filter_type == 1: # low nat = stopb / passb elif filter_type == 2: # high nat = passb / stopb elif filter_type == 3: # stop - wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0]-1e-12, - args=(0,passb,stopb,gpass,gstop,'cheby'), disp=0) + wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0] - 1e-12, + args=(0, passb, stopb, gpass, gstop, 'cheby'), + disp=0) passb[0] = wp0 - wp1 = optimize.fminbound(band_stop_obj, stopb[1]+1e-12, passb[1], - args=(1,passb,stopb,gpass,gstop,'cheby'), disp=0) + wp1 = optimize.fminbound(band_stop_obj, stopb[1] + 1e-12, passb[1], + args=(1, passb, stopb, gpass, gstop, 'cheby'), + disp=0) passb[1] = wp1 - nat = (stopb * (passb[0]-passb[1])) / (stopb**2 - passb[0]*passb[1]) + nat = ((stopb * (passb[0] - passb[1])) / + (stopb ** 2 - passb[0] * passb[1])) elif filter_type == 4: # pass - nat = (stopb**2 - passb[0]*passb[1]) / (stopb* (passb[0]-passb[1])) + nat = ((stopb ** 2 - passb[0] * passb[1]) / + (stopb * (passb[0] - passb[1]))) nat = min(abs(nat)) - GSTOP = 10**(0.1*abs(gstop)) - GPASS = 10**(0.1*abs(gpass)) - ord = int(ceil(arccosh(sqrt((GSTOP-1.0) / (GPASS-1.0))) / arccosh(nat))) + GSTOP = 10 ** (0.1 * abs(gstop)) + GPASS = 10 ** (0.1 * abs(gpass)) + ord = int(ceil(arccosh(sqrt((GSTOP - 1.0) / (GPASS - 1.0))) / + arccosh(nat))) # Natural frequencies are just the passband edges if not analog: - wn = (2.0/pi)*arctan(passb) + wn = (2.0 / pi) * arctan(passb) else: wn = passb @@ -925,8 +973,8 @@ def cheb2ord(wp, ws, gpass, gstop, analog=0): Description: Return the order of the lowest order digital Chebyshev Type II filter - that loses no more than gpass dB in the passband and has at least gstop dB - attenuation in the stopband. + that loses no more than gpass dB in the passband and has at least + gstop dB attenuation in the stopband. Parameters ---------- @@ -958,7 +1006,7 @@ def cheb2ord(wp, ws, gpass, gstop, analog=0): """ wp = atleast_1d(wp) ws = atleast_1d(ws) - filter_type = 2*(len(wp)-1) + filter_type = 2 * (len(wp) - 1) if wp[0] < ws[0]: filter_type += 1 else: @@ -966,39 +1014,42 @@ def cheb2ord(wp, ws, gpass, gstop, analog=0): # Pre-wagpass frequencies if not analog: - passb = tan(pi*wp/2.0) - stopb = tan(pi*ws/2.0) + passb = tan(pi * wp / 2.0) + stopb = tan(pi * ws / 2.0) else: - passb = wp*1.0 - stopb = ws*1.0 + passb = wp * 1.0 + stopb = ws * 1.0 if filter_type == 1: # low nat = stopb / passb elif filter_type == 2: # high nat = passb / stopb elif filter_type == 3: # stop - wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0]-1e-12, - args=(0,passb,stopb,gpass,gstop,'cheby'), + wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0] - 1e-12, + args=(0, passb, stopb, gpass, gstop, 'cheby'), disp=0) passb[0] = wp0 - wp1 = optimize.fminbound(band_stop_obj, stopb[1]+1e-12, passb[1], - args=(1,passb,stopb,gpass,gstop,'cheby'), + wp1 = optimize.fminbound(band_stop_obj, stopb[1] + 1e-12, passb[1], + args=(1, passb, stopb, gpass, gstop, 'cheby'), disp=0) passb[1] = wp1 - nat = (stopb * (passb[0]-passb[1])) / (stopb**2 - passb[0]*passb[1]) + nat = ((stopb * (passb[0] - passb[1])) / + (stopb ** 2 - passb[0] * passb[1])) elif filter_type == 4: # pass - nat = (stopb**2 - passb[0]*passb[1]) / (stopb* (passb[0]-passb[1])) + nat = ((stopb ** 2 - passb[0] * passb[1]) / + (stopb * (passb[0] - passb[1]))) nat = min(abs(nat)) - GSTOP = 10**(0.1*abs(gstop)) - GPASS = 10**(0.1*abs(gpass)) - ord = int(ceil(arccosh(sqrt((GSTOP-1.0) / (GPASS-1.0))) / arccosh(nat))) + GSTOP = 10 ** (0.1 * abs(gstop)) + GPASS = 10 ** (0.1 * abs(gpass)) + ord = int(ceil(arccosh(sqrt((GSTOP - 1.0) / (GPASS - 1.0))) / + arccosh(nat))) # Find frequency where analog response is -gpass dB. # Then convert back from low-pass prototype to the original filter. - new_freq = cosh(1.0/ord * arccosh(sqrt((GSTOP-1.0)/(GPASS-1.0)))) + new_freq = cosh(1.0 / ord * arccosh(sqrt((GSTOP - 1.0) / (GPASS - 1.0)))) new_freq = 1.0 / new_freq if filter_type == 1: @@ -1006,20 +1057,20 @@ def cheb2ord(wp, ws, gpass, gstop, analog=0): elif filter_type == 2: nat = passb * new_freq elif filter_type == 3: - nat = numpy.zeros(2,float) - nat[0] = new_freq / 2.0 * (passb[0]-passb[1]) + \ - sqrt(new_freq**2 * (passb[1]-passb[0])**2 / 4.0 + \ - passb[1] * passb[0]) + nat = numpy.zeros(2, float) + nat[0] = (new_freq / 2.0 * (passb[0] - passb[1]) + + sqrt(new_freq ** 2 * (passb[1] - passb[0]) ** 2 / 4.0 + + passb[1] * passb[0])) nat[1] = passb[1] * passb[0] / nat[0] elif filter_type == 4: - nat = numpy.zeros(2,float) - nat[0] = 1.0/(2.0*new_freq) * (passb[0] - passb[1]) + \ - sqrt((passb[1]-passb[0])**2 / (4.0*new_freq**2) + \ - passb[1] * passb[0]) + nat = numpy.zeros(2, float) + nat[0] = (1.0 / (2.0 * new_freq) * (passb[0] - passb[1]) + + sqrt((passb[1] - passb[0]) ** 2 / (4.0 * new_freq ** 2) + + passb[1] * passb[0])) nat[1] = passb[0] * passb[1] / nat[0] if not analog: - wn = (2.0/pi)*arctan(nat) + wn = (2.0 / pi) * arctan(nat) else: wn = nat @@ -1065,48 +1116,50 @@ def ellipord(wp, ws, gpass, gstop, analog=0): """ wp = atleast_1d(wp) ws = atleast_1d(ws) - filter_type = 2*(len(wp)-1) + filter_type = 2 * (len(wp) - 1) filter_type += 1 if wp[0] >= ws[0]: filter_type += 1 # Pre-wagpass frequencies if analog: - passb = wp*1.0 - stopb = ws*1.0 + passb = wp * 1.0 + stopb = ws * 1.0 else: - passb = tan(wp*pi/2.0) - stopb = tan(ws*pi/2.0) + passb = tan(wp * pi / 2.0) + stopb = tan(ws * pi / 2.0) if filter_type == 1: # low nat = stopb / passb elif filter_type == 2: # high nat = passb / stopb elif filter_type == 3: # stop - wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0]-1e-12, - args=(0,passb,stopb,gpass,gstop,'ellip'), + wp0 = optimize.fminbound(band_stop_obj, passb[0], stopb[0] - 1e-12, + args=(0, passb, stopb, gpass, gstop, 'ellip'), disp=0) passb[0] = wp0 - wp1 = optimize.fminbound(band_stop_obj, stopb[1]+1e-12, passb[1], - args=(1,passb,stopb,gpass,gstop,'ellip'), + wp1 = optimize.fminbound(band_stop_obj, stopb[1] + 1e-12, passb[1], + args=(1, passb, stopb, gpass, gstop, 'ellip'), disp=0) passb[1] = wp1 - nat = (stopb * (passb[0]-passb[1])) / (stopb**2 - passb[0]*passb[1]) + nat = ((stopb * (passb[0] - passb[1])) / + (stopb ** 2 - passb[0] * passb[1])) elif filter_type == 4: # pass - nat = (stopb**2 - passb[0]*passb[1]) / (stopb* (passb[0]-passb[1])) + nat = ((stopb ** 2 - passb[0] * passb[1]) / + (stopb * (passb[0] - passb[1]))) nat = min(abs(nat)) - GSTOP = 10**(0.1*gstop) - GPASS = 10**(0.1*gpass) - arg1 = sqrt( (GPASS-1.0) / (GSTOP-1.0) ) + GSTOP = 10 ** (0.1 * gstop) + GPASS = 10 ** (0.1 * gpass) + arg1 = sqrt((GPASS - 1.0) / (GSTOP - 1.0)) arg0 = 1.0 / nat - d0 = special.ellipk([arg0**2, 1-arg0**2]) - d1 = special.ellipk([arg1**2, 1-arg1**2]) - ord = int(ceil(d0[0]*d1[1] / (d0[1]*d1[0]))) + d0 = special.ellipk([arg0 ** 2, 1 - arg0 ** 2]) + d1 = special.ellipk([arg1 ** 2, 1 - arg1 ** 2]) + ord = int(ceil(d0[0] * d1[1] / (d0[1] * d1[0]))) if not analog: - wn = arctan(passb)*2.0/pi + wn = arctan(passb) * 2.0 / pi else: wn = passb @@ -1114,60 +1167,66 @@ def ellipord(wp, ws, gpass, gstop, analog=0): wn = wn[0] return ord, wn + def buttap(N): """Return (z,p,k) zero, pole, gain for analog prototype of an Nth order Butterworth filter.""" z = [] - n = numpy.arange(1,N+1) - p = numpy.exp(1j*(2*n-1)/(2.0*N)*pi)*1j + n = numpy.arange(1, N + 1) + p = numpy.exp(1j * (2 * n - 1) / (2.0 * N) * pi) * 1j k = 1 return z, p, k + def cheb1ap(N, rp): """Return (z,p,k) zero, pole, gain for Nth order Chebyshev type I lowpass analog filter prototype with `rp` decibels of ripple in the passband. """ z = [] - eps = numpy.sqrt(10**(0.1*rp)-1.0) - n = numpy.arange(1,N+1) - mu = 1.0/N * numpy.log((1.0+numpy.sqrt(1+eps*eps)) / eps) - theta = pi/2.0 * (2*n-1.0)/N - p = -numpy.sinh(mu)*numpy.sin(theta) + 1j*numpy.cosh(mu)*numpy.cos(theta) - k = numpy.prod(-p,axis=0).real + eps = numpy.sqrt(10 ** (0.1 * rp) - 1.0) + n = numpy.arange(1, N + 1) + mu = 1.0 / N * numpy.log((1.0 + numpy.sqrt(1 + eps * eps)) / eps) + theta = pi / 2.0 * (2 * n - 1.0) / N + p = (-numpy.sinh(mu) * numpy.sin(theta) + + 1j * numpy.cosh(mu) * numpy.cos(theta)) + k = numpy.prod(-p, axis=0).real if N % 2 == 0: - k = k / sqrt((1+eps*eps)) + k = k / sqrt((1 + eps * eps)) return z, p, k - pass + def cheb2ap(N, rs): """Return (z,p,k) zero, pole, gain for Nth order Chebyshev type II lowpass analog filter prototype with `rs` decibels of ripple in the stopband. """ - de = 1.0/sqrt(10**(0.1*rs)-1) - mu = arcsinh(1.0/de)/N + de = 1.0 / sqrt(10 ** (0.1 * rs) - 1) + mu = arcsinh(1.0 / de) / N if N % 2: m = N - 1 - n = numpy.concatenate((numpy.arange(1,N-1,2),numpy.arange(N+2,2*N,2))) + n = numpy.concatenate((numpy.arange(1, N - 1, 2), + numpy.arange(N + 2, 2 * N, 2))) else: m = N - n = numpy.arange(1,2*N,2) + n = numpy.arange(1, 2 * N, 2) - z = conjugate(1j / cos(n*pi/(2.0*N))) - p = exp(1j*(pi*numpy.arange(1,2*N,2)/(2.0*N) + pi/2.0)) - p = sinh(mu) * p.real + 1j*cosh(mu)*p.imag + z = conjugate(1j / cos(n * pi / (2.0 * N))) + p = exp(1j * (pi * numpy.arange(1, 2 * N, 2) / (2.0 * N) + pi / 2.0)) + p = sinh(mu) * p.real + 1j * cosh(mu) * p.imag p = 1.0 / p - k = (numpy.prod(-p,axis=0)/numpy.prod(-z,axis=0)).real + k = (numpy.prod(-p, axis=0) / numpy.prod(-z, axis=0)).real return z, p, k EPSILON = 2e-16 + def _vratio(u, ineps, mp): - [s,c,d,phi] = special.ellipj(u,mp) - ret = abs(ineps - s/c) + [s, c, d, phi] = special.ellipj(u, mp) + ret = abs(ineps - s / c) return ret + def _kratio(m, k_ratio): m = float(m) if m < 0: @@ -1175,7 +1234,7 @@ def _kratio(m, k_ratio): if m > 1: m = 1.0 if abs(m) > EPSILON and (abs(m) + EPSILON) < 1: - k = special.ellipk([m,1-m]) + k = special.ellipk([m, 1 - m]) r = k[0] / k[1] - k_ratio elif abs(m) > EPSILON: r = -k_ratio @@ -1183,6 +1242,7 @@ def _kratio(m, k_ratio): r = 1e20 return abs(r) + def ellipap(N, rp, rs): """Return (z,p,k) zeros, poles, and gain of an Nth order normalized prototype elliptic analog lowpass filter with `rp` decibels of ripple in @@ -1195,23 +1255,24 @@ def ellipap(N, rp, rs): """ if N == 1: - p = -sqrt(1.0/(10**(0.1*rp)-1.0)) + p = -sqrt(1.0 / (10 ** (0.1 * rp) - 1.0)) k = -p z = [] return z, p, k - eps = numpy.sqrt(10**(0.1*rp)-1) - ck1 = eps / numpy.sqrt(10**(0.1*rs)-1) - ck1p = numpy.sqrt(1-ck1*ck1) + eps = numpy.sqrt(10 ** (0.1 * rp) - 1) + ck1 = eps / numpy.sqrt(10 ** (0.1 * rs) - 1) + ck1p = numpy.sqrt(1 - ck1 * ck1) if ck1p == 1: - raise ValueError("Cannot design a filter with given rp and rs specifications.") + raise ValueError("Cannot design a filter with given rp and rs" + " specifications.") wp = 1 - val = special.ellipk([ck1*ck1,ck1p*ck1p]) - if abs(1-ck1p*ck1p) < EPSILON: + val = special.ellipk([ck1 * ck1, ck1p * ck1p]) + if abs(1 - ck1p * ck1p) < EPSILON: krat = 0 else: - krat = N*val[0] / val[1] + krat = N * val[0] / val[1] m = optimize.fmin(_kratio, [0.5], args=(krat,), maxfun=250, maxiter=250, disp=0) @@ -1221,33 +1282,36 @@ def ellipap(N, rp, rs): capk = special.ellipk(m) ws = wp / sqrt(m) - m1 = 1-m + m1 = 1 - m - j = numpy.arange(1-N%2,N,2) + j = numpy.arange(1 - N % 2, N, 2) jj = len(j) - [s,c,d,phi] = special.ellipj(j*capk/N,m*numpy.ones(jj)) - snew = numpy.compress(abs(s) > EPSILON, s,axis=-1) - z = 1.0 / (sqrt(m)*snew) - z = 1j*z - z = numpy.concatenate((z,conjugate(z))) + [s, c, d, phi] = special.ellipj(j * capk / N, m * numpy.ones(jj)) + snew = numpy.compress(abs(s) > EPSILON, s, axis=-1) + z = 1.0 / (sqrt(m) * snew) + z = 1j * z + z = numpy.concatenate((z, conjugate(z))) - r = optimize.fmin(_vratio, special.ellipk(m), args=(1./eps, ck1p*ck1p), + r = optimize.fmin(_vratio, special.ellipk(m), args=(1. / eps, ck1p * ck1p), maxfun=250, maxiter=250, disp=0) - v0 = capk * r / (N*val[0]) + v0 = capk * r / (N * val[0]) - [sv,cv,dv,phi] = special.ellipj(v0,1-m) - p = -(c*d*sv*cv + 1j*s*dv) / (1-(d*sv)**2.0) + [sv, cv, dv, phi] = special.ellipj(v0, 1 - m) + p = -(c * d * sv * cv + 1j * s * dv) / (1 - (d * sv) ** 2.0) if N % 2: - newp = numpy.compress(abs(p.imag) > EPSILON*numpy.sqrt(numpy.sum(p*numpy.conjugate(p),axis=0).real), p,axis=-1) - p = numpy.concatenate((p,conjugate(newp))) + newp = numpy.compress(abs(p.imag) > EPSILON * + numpy.sqrt(numpy.sum(p * numpy.conjugate(p), + axis=0).real), + p, axis=-1) + p = numpy.concatenate((p, conjugate(newp))) else: - p = numpy.concatenate((p,conjugate(p))) + p = numpy.concatenate((p, conjugate(p))) - k = (numpy.prod(-p,axis=0) / numpy.prod(-z,axis=0)).real + k = (numpy.prod(-p, axis=0) / numpy.prod(-z, axis=0)).real if N % 2 == 0: - k = k / numpy.sqrt((1+eps*eps)) + k = k / numpy.sqrt((1 + eps * eps)) return z, p, k @@ -1258,390 +1322,390 @@ def besselap(N): z = [] k = 1 if N == 0: - p = []; + p = [] elif N == 1: p = [-1] elif N == 2: - p = [-.8660254037844386467637229+.4999999999999999999999996*1j, - -.8660254037844386467637229-.4999999999999999999999996*1j] + p = [-.8660254037844386467637229 + .4999999999999999999999996j, + -.8660254037844386467637229 - .4999999999999999999999996j] elif N == 3: p = [-.9416000265332067855971980, - -.7456403858480766441810907-.7113666249728352680992154*1j, - -.7456403858480766441810907+.7113666249728352680992154*1j] + -.7456403858480766441810907 - .7113666249728352680992154j, + -.7456403858480766441810907 + .7113666249728352680992154j] elif N == 4: - p = [-.6572111716718829545787781-.8301614350048733772399715*1j, - -.6572111716718829545787788+.8301614350048733772399715*1j, - -.9047587967882449459642637-.2709187330038746636700923*1j, - -.9047587967882449459642624+.2709187330038746636700926*1j] + p = [-.6572111716718829545787781 - .8301614350048733772399715j, + -.6572111716718829545787788 + .8301614350048733772399715j, + -.9047587967882449459642637 - .2709187330038746636700923j, + -.9047587967882449459642624 + .2709187330038746636700926j] elif N == 5: p = [-.9264420773877602247196260, - -.8515536193688395541722677-.4427174639443327209850002*1j, - -.8515536193688395541722677+.4427174639443327209850002*1j, - -.5905759446119191779319432-.9072067564574549539291747*1j, - -.5905759446119191779319432+.9072067564574549539291747*1j] + -.8515536193688395541722677 - .4427174639443327209850002j, + -.8515536193688395541722677 + .4427174639443327209850002j, + -.5905759446119191779319432 - .9072067564574549539291747j, + -.5905759446119191779319432 + .9072067564574549539291747j] elif N == 6: - p = [-.9093906830472271808050953-.1856964396793046769246397*1j, - -.9093906830472271808050953+.1856964396793046769246397*1j, - -.7996541858328288520243325-.5621717346937317988594118*1j, - -.7996541858328288520243325+.5621717346937317988594118*1j, - -.5385526816693109683073792-.9616876881954277199245657*1j, - -.5385526816693109683073792+.9616876881954277199245657*1j] + p = [-.9093906830472271808050953 - .1856964396793046769246397j, + -.9093906830472271808050953 + .1856964396793046769246397j, + -.7996541858328288520243325 - .5621717346937317988594118j, + -.7996541858328288520243325 + .5621717346937317988594118j, + -.5385526816693109683073792 - .9616876881954277199245657j, + -.5385526816693109683073792 + .9616876881954277199245657j] elif N == 7: p = [-.9194871556490290014311619, - -.8800029341523374639772340-.3216652762307739398381830*1j, - -.8800029341523374639772340+.3216652762307739398381830*1j, - -.7527355434093214462291616-.6504696305522550699212995*1j, - -.7527355434093214462291616+.6504696305522550699212995*1j, - -.4966917256672316755024763-1.002508508454420401230220*1j, - -.4966917256672316755024763+1.002508508454420401230220*1j] + -.8800029341523374639772340 - .3216652762307739398381830j, + -.8800029341523374639772340 + .3216652762307739398381830j, + -.7527355434093214462291616 - .6504696305522550699212995j, + -.7527355434093214462291616 + .6504696305522550699212995j, + -.4966917256672316755024763 - 1.002508508454420401230220j, + -.4966917256672316755024763 + 1.002508508454420401230220j] elif N == 8: - p = [-.9096831546652910216327629-.1412437976671422927888150*1j, - -.9096831546652910216327629+.1412437976671422927888150*1j, - -.8473250802359334320103023-.4259017538272934994996429*1j, - -.8473250802359334320103023+.4259017538272934994996429*1j, - -.7111381808485399250796172-.7186517314108401705762571*1j, - -.7111381808485399250796172+.7186517314108401705762571*1j, - -.4621740412532122027072175-1.034388681126901058116589*1j, - -.4621740412532122027072175+1.034388681126901058116589*1j] + p = [-.9096831546652910216327629 - .1412437976671422927888150j, + -.9096831546652910216327629 + .1412437976671422927888150j, + -.8473250802359334320103023 - .4259017538272934994996429j, + -.8473250802359334320103023 + .4259017538272934994996429j, + -.7111381808485399250796172 - .7186517314108401705762571j, + -.7111381808485399250796172 + .7186517314108401705762571j, + -.4621740412532122027072175 - 1.034388681126901058116589j, + -.4621740412532122027072175 + 1.034388681126901058116589j] elif N == 9: p = [-.9154957797499037686769223, - -.8911217017079759323183848-.2526580934582164192308115*1j, - -.8911217017079759323183848+.2526580934582164192308115*1j, - -.8148021112269012975514135-.5085815689631499483745341*1j, - -.8148021112269012975514135+.5085815689631499483745341*1j, - -.6743622686854761980403401-.7730546212691183706919682*1j, - -.6743622686854761980403401+.7730546212691183706919682*1j, - -.4331415561553618854685942-1.060073670135929666774323*1j, - -.4331415561553618854685942+1.060073670135929666774323*1j] + -.8911217017079759323183848 - .2526580934582164192308115j, + -.8911217017079759323183848 + .2526580934582164192308115j, + -.8148021112269012975514135 - .5085815689631499483745341j, + -.8148021112269012975514135 + .5085815689631499483745341j, + -.6743622686854761980403401 - .7730546212691183706919682j, + -.6743622686854761980403401 + .7730546212691183706919682j, + -.4331415561553618854685942 - 1.060073670135929666774323j, + -.4331415561553618854685942 + 1.060073670135929666774323j] elif N == 10: - p = [-.9091347320900502436826431-.1139583137335511169927714*1j, - -.9091347320900502436826431+.1139583137335511169927714*1j, - -.8688459641284764527921864-.3430008233766309973110589*1j, - -.8688459641284764527921864+.3430008233766309973110589*1j, - -.7837694413101441082655890-.5759147538499947070009852*1j, - -.7837694413101441082655890+.5759147538499947070009852*1j, - -.6417513866988316136190854-.8175836167191017226233947*1j, - -.6417513866988316136190854+.8175836167191017226233947*1j, - -.4083220732868861566219785-1.081274842819124562037210*1j, - -.4083220732868861566219785+1.081274842819124562037210*1j] + p = [-.9091347320900502436826431 - .1139583137335511169927714j, + -.9091347320900502436826431 + .1139583137335511169927714j, + -.8688459641284764527921864 - .3430008233766309973110589j, + -.8688459641284764527921864 + .3430008233766309973110589j, + -.7837694413101441082655890 - .5759147538499947070009852j, + -.7837694413101441082655890 + .5759147538499947070009852j, + -.6417513866988316136190854 - .8175836167191017226233947j, + -.6417513866988316136190854 + .8175836167191017226233947j, + -.4083220732868861566219785 - 1.081274842819124562037210j, + -.4083220732868861566219785 + 1.081274842819124562037210j] elif N == 11: p = [-.9129067244518981934637318, - -.8963656705721166099815744-.2080480375071031919692341*1j - -.8963656705721166099815744+.2080480375071031919692341*1j, - -.8453044014712962954184557-.4178696917801248292797448*1j, - -.8453044014712962954184557+.4178696917801248292797448*1j, - -.7546938934722303128102142-.6319150050721846494520941*1j, - -.7546938934722303128102142+.6319150050721846494520941*1j, - -.6126871554915194054182909-.8547813893314764631518509*1j, - -.6126871554915194054182909+.8547813893314764631518509*1j, - -.3868149510055090879155425-1.099117466763120928733632*1j, - -.3868149510055090879155425+1.099117466763120928733632*1j] + -.8963656705721166099815744 - .2080480375071031919692341j, + -.8963656705721166099815744 + .2080480375071031919692341j, + -.8453044014712962954184557 - .4178696917801248292797448j, + -.8453044014712962954184557 + .4178696917801248292797448j, + -.7546938934722303128102142 - .6319150050721846494520941j, + -.7546938934722303128102142 + .6319150050721846494520941j, + -.6126871554915194054182909 - .8547813893314764631518509j, + -.6126871554915194054182909 + .8547813893314764631518509j, + -.3868149510055090879155425 - 1.099117466763120928733632j, + -.3868149510055090879155425 + 1.099117466763120928733632j] elif N == 12: - p = [-.9084478234140682638817772-95506365213450398415258360.0e-27*1j, - -.9084478234140682638817772+95506365213450398415258360.0e-27*1j, - -.8802534342016826507901575-.2871779503524226723615457*1j, - -.8802534342016826507901575+.2871779503524226723615457*1j, - -.8217296939939077285792834-.4810212115100676440620548*1j, - -.8217296939939077285792834+.4810212115100676440620548*1j, - -.7276681615395159454547013-.6792961178764694160048987*1j, - -.7276681615395159454547013+.6792961178764694160048987*1j, - -.5866369321861477207528215-.8863772751320727026622149*1j, - -.5866369321861477207528215+.8863772751320727026622149*1j, - -.3679640085526312839425808-1.114373575641546257595657*1j, - -.3679640085526312839425808+1.114373575641546257595657*1j] + p = [-.9084478234140682638817772 - 95506365213450398415258360.0e-27j, + -.9084478234140682638817772 + 95506365213450398415258360.0e-27j, + -.8802534342016826507901575 - .2871779503524226723615457j, + -.8802534342016826507901575 + .2871779503524226723615457j, + -.8217296939939077285792834 - .4810212115100676440620548j, + -.8217296939939077285792834 + .4810212115100676440620548j, + -.7276681615395159454547013 - .6792961178764694160048987j, + -.7276681615395159454547013 + .6792961178764694160048987j, + -.5866369321861477207528215 - .8863772751320727026622149j, + -.5866369321861477207528215 + .8863772751320727026622149j, + -.3679640085526312839425808 - 1.114373575641546257595657j, + -.3679640085526312839425808 + 1.114373575641546257595657j] elif N == 13: p = [-.9110914665984182781070663, - -.8991314665475196220910718-.1768342956161043620980863*1j, - -.8991314665475196220910718+.1768342956161043620980863*1j, - -.8625094198260548711573628-.3547413731172988997754038*1j, - -.8625094198260548711573628+.3547413731172988997754038*1j, - -.7987460692470972510394686-.5350752120696801938272504*1j, - -.7987460692470972510394686+.5350752120696801938272504*1j, - -.7026234675721275653944062-.7199611890171304131266374*1j, - -.7026234675721275653944062+.7199611890171304131266374*1j, - -.5631559842430199266325818-.9135900338325109684927731*1j, - -.5631559842430199266325818+.9135900338325109684927731*1j, - -.3512792323389821669401925-1.127591548317705678613239*1j, - -.3512792323389821669401925+1.127591548317705678613239*1j] + -.8991314665475196220910718 - .1768342956161043620980863j, + -.8991314665475196220910718 + .1768342956161043620980863j, + -.8625094198260548711573628 - .3547413731172988997754038j, + -.8625094198260548711573628 + .3547413731172988997754038j, + -.7987460692470972510394686 - .5350752120696801938272504j, + -.7987460692470972510394686 + .5350752120696801938272504j, + -.7026234675721275653944062 - .7199611890171304131266374j, + -.7026234675721275653944062 + .7199611890171304131266374j, + -.5631559842430199266325818 - .9135900338325109684927731j, + -.5631559842430199266325818 + .9135900338325109684927731j, + -.3512792323389821669401925 - 1.127591548317705678613239j, + -.3512792323389821669401925 + 1.127591548317705678613239j] elif N == 14: - p = [-.9077932138396487614720659-82196399419401501888968130.0e-27*1j, - -.9077932138396487614720659+82196399419401501888968130.0e-27*1j, - -.8869506674916445312089167-.2470079178765333183201435*1j, - -.8869506674916445312089167+.2470079178765333183201435*1j, - -.8441199160909851197897667-.4131653825102692595237260*1j, - -.8441199160909851197897667+.4131653825102692595237260*1j, - -.7766591387063623897344648-.5819170677377608590492434*1j, - -.7766591387063623897344648+.5819170677377608590492434*1j, - -.6794256425119233117869491-.7552857305042033418417492*1j, - -.6794256425119233117869491+.7552857305042033418417492*1j, - -.5418766775112297376541293-.9373043683516919569183099*1j, - -.5418766775112297376541293+.9373043683516919569183099*1j, - -.3363868224902037330610040-1.139172297839859991370924*1j, - -.3363868224902037330610040+1.139172297839859991370924*1j] + p = [-.9077932138396487614720659 - 82196399419401501888968130.0e-27j, + -.9077932138396487614720659 + 82196399419401501888968130.0e-27j, + -.8869506674916445312089167 - .2470079178765333183201435j, + -.8869506674916445312089167 + .2470079178765333183201435j, + -.8441199160909851197897667 - .4131653825102692595237260j, + -.8441199160909851197897667 + .4131653825102692595237260j, + -.7766591387063623897344648 - .5819170677377608590492434j, + -.7766591387063623897344648 + .5819170677377608590492434j, + -.6794256425119233117869491 - .7552857305042033418417492j, + -.6794256425119233117869491 + .7552857305042033418417492j, + -.5418766775112297376541293 - .9373043683516919569183099j, + -.5418766775112297376541293 + .9373043683516919569183099j, + -.3363868224902037330610040 - 1.139172297839859991370924j, + -.3363868224902037330610040 + 1.139172297839859991370924j] elif N == 15: p = [-.9097482363849064167228581, - -.9006981694176978324932918-.1537681197278439351298882*1j, - -.9006981694176978324932918+.1537681197278439351298882*1j, - -.8731264620834984978337843-.3082352470564267657715883*1j, - -.8731264620834984978337843+.3082352470564267657715883*1j, - -.8256631452587146506294553-.4642348752734325631275134*1j, - -.8256631452587146506294553+.4642348752734325631275134*1j, - -.7556027168970728127850416-.6229396358758267198938604*1j, - -.7556027168970728127850416+.6229396358758267198938604*1j, - -.6579196593110998676999362-.7862895503722515897065645*1j, - -.6579196593110998676999362+.7862895503722515897065645*1j, - -.5224954069658330616875186-.9581787261092526478889345*1j, - -.5224954069658330616875186+.9581787261092526478889345*1j, - -.3229963059766444287113517-1.149416154583629539665297*1j, - -.3229963059766444287113517+1.149416154583629539665297*1j] + -.9006981694176978324932918 - .1537681197278439351298882j, + -.9006981694176978324932918 + .1537681197278439351298882j, + -.8731264620834984978337843 - .3082352470564267657715883j, + -.8731264620834984978337843 + .3082352470564267657715883j, + -.8256631452587146506294553 - .4642348752734325631275134j, + -.8256631452587146506294553 + .4642348752734325631275134j, + -.7556027168970728127850416 - .6229396358758267198938604j, + -.7556027168970728127850416 + .6229396358758267198938604j, + -.6579196593110998676999362 - .7862895503722515897065645j, + -.6579196593110998676999362 + .7862895503722515897065645j, + -.5224954069658330616875186 - .9581787261092526478889345j, + -.5224954069658330616875186 + .9581787261092526478889345j, + -.3229963059766444287113517 - 1.149416154583629539665297j, + -.3229963059766444287113517 + 1.149416154583629539665297j] elif N == 16: - p = [-.9072099595087001356491337-72142113041117326028823950.0e-27*1j, - -.9072099595087001356491337+72142113041117326028823950.0e-27*1j, - -.8911723070323647674780132-.2167089659900576449410059*1j, - -.8911723070323647674780132+.2167089659900576449410059*1j, - -.8584264231521330481755780-.3621697271802065647661080*1j, - -.8584264231521330481755780+.3621697271802065647661080*1j, - -.8074790293236003885306146-.5092933751171800179676218*1j, - -.8074790293236003885306146+.5092933751171800179676218*1j, - -.7356166304713115980927279-.6591950877860393745845254*1j, - -.7356166304713115980927279+.6591950877860393745845254*1j, - -.6379502514039066715773828-.8137453537108761895522580*1j, - -.6379502514039066715773828+.8137453537108761895522580*1j, - -.5047606444424766743309967-.9767137477799090692947061*1j, - -.5047606444424766743309967+.9767137477799090692947061*1j, - -.3108782755645387813283867-1.158552841199330479412225*1j, - -.3108782755645387813283867+1.158552841199330479412225*1j] + p = [-.9072099595087001356491337 - 72142113041117326028823950.0e-27j, + -.9072099595087001356491337 + 72142113041117326028823950.0e-27j, + -.8911723070323647674780132 - .2167089659900576449410059j, + -.8911723070323647674780132 + .2167089659900576449410059j, + -.8584264231521330481755780 - .3621697271802065647661080j, + -.8584264231521330481755780 + .3621697271802065647661080j, + -.8074790293236003885306146 - .5092933751171800179676218j, + -.8074790293236003885306146 + .5092933751171800179676218j, + -.7356166304713115980927279 - .6591950877860393745845254j, + -.7356166304713115980927279 + .6591950877860393745845254j, + -.6379502514039066715773828 - .8137453537108761895522580j, + -.6379502514039066715773828 + .8137453537108761895522580j, + -.5047606444424766743309967 - .9767137477799090692947061j, + -.5047606444424766743309967 + .9767137477799090692947061j, + -.3108782755645387813283867 - 1.158552841199330479412225j, + -.3108782755645387813283867 + 1.158552841199330479412225j] elif N == 17: p = [-.9087141161336397432860029, - -.9016273850787285964692844-.1360267995173024591237303*1j, - -.9016273850787285964692844+.1360267995173024591237303*1j, - -.8801100704438627158492165-.2725347156478803885651973*1j, - -.8801100704438627158492165+.2725347156478803885651973*1j, - -.8433414495836129204455491-.4100759282910021624185986*1j, - -.8433414495836129204455491+.4100759282910021624185986*1j, - -.7897644147799708220288138-.5493724405281088674296232*1j, - -.7897644147799708220288138+.5493724405281088674296232*1j, - -.7166893842372349049842743-.6914936286393609433305754*1j, - -.7166893842372349049842743+.6914936286393609433305754*1j, - -.6193710717342144521602448-.8382497252826992979368621*1j, - -.6193710717342144521602448+.8382497252826992979368621*1j, - -.4884629337672704194973683-.9932971956316781632345466*1j, - -.4884629337672704194973683+.9932971956316781632345466*1j, - -.2998489459990082015466971-1.166761272925668786676672*1j, - -.2998489459990082015466971+1.166761272925668786676672*1j] + -.9016273850787285964692844 - .1360267995173024591237303j, + -.9016273850787285964692844 + .1360267995173024591237303j, + -.8801100704438627158492165 - .2725347156478803885651973j, + -.8801100704438627158492165 + .2725347156478803885651973j, + -.8433414495836129204455491 - .4100759282910021624185986j, + -.8433414495836129204455491 + .4100759282910021624185986j, + -.7897644147799708220288138 - .5493724405281088674296232j, + -.7897644147799708220288138 + .5493724405281088674296232j, + -.7166893842372349049842743 - .6914936286393609433305754j, + -.7166893842372349049842743 + .6914936286393609433305754j, + -.6193710717342144521602448 - .8382497252826992979368621j, + -.6193710717342144521602448 + .8382497252826992979368621j, + -.4884629337672704194973683 - .9932971956316781632345466j, + -.4884629337672704194973683 + .9932971956316781632345466j, + -.2998489459990082015466971 - 1.166761272925668786676672j, + -.2998489459990082015466971 + 1.166761272925668786676672j] elif N == 18: - p = [-.9067004324162775554189031-64279241063930693839360680.0e-27*1j, - -.9067004324162775554189031+64279241063930693839360680.0e-27*1j, - -.8939764278132455733032155-.1930374640894758606940586*1j, - -.8939764278132455733032155+.1930374640894758606940586*1j, - -.8681095503628830078317207-.3224204925163257604931634*1j, - -.8681095503628830078317207+.3224204925163257604931634*1j, - -.8281885016242836608829018-.4529385697815916950149364*1j, - -.8281885016242836608829018+.4529385697815916950149364*1j, - -.7726285030739558780127746-.5852778162086640620016316*1j, - -.7726285030739558780127746+.5852778162086640620016316*1j, - -.6987821445005273020051878-.7204696509726630531663123*1j, - -.6987821445005273020051878+.7204696509726630531663123*1j, - -.6020482668090644386627299-.8602708961893664447167418*1j, - -.6020482668090644386627299+.8602708961893664447167418*1j, - -.4734268069916151511140032-1.008234300314801077034158*1j, - -.4734268069916151511140032+1.008234300314801077034158*1j, - -.2897592029880489845789953-1.174183010600059128532230*1j, - -.2897592029880489845789953+1.174183010600059128532230*1j] + p = [-.9067004324162775554189031 - 64279241063930693839360680.0e-27j, + -.9067004324162775554189031 + 64279241063930693839360680.0e-27j, + -.8939764278132455733032155 - .1930374640894758606940586j, + -.8939764278132455733032155 + .1930374640894758606940586j, + -.8681095503628830078317207 - .3224204925163257604931634j, + -.8681095503628830078317207 + .3224204925163257604931634j, + -.8281885016242836608829018 - .4529385697815916950149364j, + -.8281885016242836608829018 + .4529385697815916950149364j, + -.7726285030739558780127746 - .5852778162086640620016316j, + -.7726285030739558780127746 + .5852778162086640620016316j, + -.6987821445005273020051878 - .7204696509726630531663123j, + -.6987821445005273020051878 + .7204696509726630531663123j, + -.6020482668090644386627299 - .8602708961893664447167418j, + -.6020482668090644386627299 + .8602708961893664447167418j, + -.4734268069916151511140032 - 1.008234300314801077034158j, + -.4734268069916151511140032 + 1.008234300314801077034158j, + -.2897592029880489845789953 - 1.174183010600059128532230j, + -.2897592029880489845789953 + 1.174183010600059128532230j] elif N == 19: p = [-.9078934217899404528985092, - -.9021937639390660668922536-.1219568381872026517578164*1j, - -.9021937639390660668922536+.1219568381872026517578164*1j, - -.8849290585034385274001112-.2442590757549818229026280*1j, - -.8849290585034385274001112+.2442590757549818229026280*1j, - -.8555768765618421591093993-.3672925896399872304734923*1j, - -.8555768765618421591093993+.3672925896399872304734923*1j, - -.8131725551578197705476160-.4915365035562459055630005*1j, - -.8131725551578197705476160+.4915365035562459055630005*1j, - -.7561260971541629355231897-.6176483917970178919174173*1j, - -.7561260971541629355231897+.6176483917970178919174173*1j, - -.6818424412912442033411634-.7466272357947761283262338*1j, - -.6818424412912442033411634+.7466272357947761283262338*1j, - -.5858613321217832644813602-.8801817131014566284786759*1j, - -.5858613321217832644813602+.8801817131014566284786759*1j, - -.4595043449730988600785456-1.021768776912671221830298*1j, - -.4595043449730988600785456+1.021768776912671221830298*1j, - -.2804866851439370027628724-1.180931628453291873626003*1j, - -.2804866851439370027628724+1.180931628453291873626003*1j] + -.9021937639390660668922536 - .1219568381872026517578164j, + -.9021937639390660668922536 + .1219568381872026517578164j, + -.8849290585034385274001112 - .2442590757549818229026280j, + -.8849290585034385274001112 + .2442590757549818229026280j, + -.8555768765618421591093993 - .3672925896399872304734923j, + -.8555768765618421591093993 + .3672925896399872304734923j, + -.8131725551578197705476160 - .4915365035562459055630005j, + -.8131725551578197705476160 + .4915365035562459055630005j, + -.7561260971541629355231897 - .6176483917970178919174173j, + -.7561260971541629355231897 + .6176483917970178919174173j, + -.6818424412912442033411634 - .7466272357947761283262338j, + -.6818424412912442033411634 + .7466272357947761283262338j, + -.5858613321217832644813602 - .8801817131014566284786759j, + -.5858613321217832644813602 + .8801817131014566284786759j, + -.4595043449730988600785456 - 1.021768776912671221830298j, + -.4595043449730988600785456 + 1.021768776912671221830298j, + -.2804866851439370027628724 - 1.180931628453291873626003j, + -.2804866851439370027628724 + 1.180931628453291873626003j] elif N == 20: - p = [-.9062570115576771146523497-57961780277849516990208850.0e-27*1j, - -.9062570115576771146523497+57961780277849516990208850.0e-27*1j, - -.8959150941925768608568248-.1740317175918705058595844*1j, - -.8959150941925768608568248+.1740317175918705058595844*1j, - -.8749560316673332850673214-.2905559296567908031706902*1j, - -.8749560316673332850673214+.2905559296567908031706902*1j, - -.8427907479956670633544106-.4078917326291934082132821*1j, - -.8427907479956670633544106+.4078917326291934082132821*1j, - -.7984251191290606875799876-.5264942388817132427317659*1j, - -.7984251191290606875799876+.5264942388817132427317659*1j, - -.7402780309646768991232610-.6469975237605228320268752*1j, - -.7402780309646768991232610+.6469975237605228320268752*1j, - -.6658120544829934193890626-.7703721701100763015154510*1j, - -.6658120544829934193890626+.7703721701100763015154510*1j, - -.5707026806915714094398061-.8982829066468255593407161*1j, - -.5707026806915714094398061+.8982829066468255593407161*1j, - -.4465700698205149555701841-1.034097702560842962315411*1j, - -.4465700698205149555701841+1.034097702560842962315411*1j, - -.2719299580251652601727704-1.187099379810885886139638*1j, - -.2719299580251652601727704+1.187099379810885886139638*1j] + p = [-.9062570115576771146523497 - 57961780277849516990208850.0e-27j, + -.9062570115576771146523497 + 57961780277849516990208850.0e-27j, + -.8959150941925768608568248 - .1740317175918705058595844j, + -.8959150941925768608568248 + .1740317175918705058595844j, + -.8749560316673332850673214 - .2905559296567908031706902j, + -.8749560316673332850673214 + .2905559296567908031706902j, + -.8427907479956670633544106 - .4078917326291934082132821j, + -.8427907479956670633544106 + .4078917326291934082132821j, + -.7984251191290606875799876 - .5264942388817132427317659j, + -.7984251191290606875799876 + .5264942388817132427317659j, + -.7402780309646768991232610 - .6469975237605228320268752j, + -.7402780309646768991232610 + .6469975237605228320268752j, + -.6658120544829934193890626 - .7703721701100763015154510j, + -.6658120544829934193890626 + .7703721701100763015154510j, + -.5707026806915714094398061 - .8982829066468255593407161j, + -.5707026806915714094398061 + .8982829066468255593407161j, + -.4465700698205149555701841 - 1.034097702560842962315411j, + -.4465700698205149555701841 + 1.034097702560842962315411j, + -.2719299580251652601727704 - 1.187099379810885886139638j, + -.2719299580251652601727704 + 1.187099379810885886139638j] elif N == 21: p = [-.9072262653142957028884077, - -.9025428073192696303995083-.1105252572789856480992275*1j, - -.9025428073192696303995083+.1105252572789856480992275*1j, - -.8883808106664449854431605-.2213069215084350419975358*1j, - -.8883808106664449854431605+.2213069215084350419975358*1j, - -.8643915813643204553970169-.3326258512522187083009453*1j, - -.8643915813643204553970169+.3326258512522187083009453*1j, - -.8299435470674444100273463-.4448177739407956609694059*1j, - -.8299435470674444100273463+.4448177739407956609694059*1j, - -.7840287980408341576100581-.5583186348022854707564856*1j, - -.7840287980408341576100581+.5583186348022854707564856*1j, - -.7250839687106612822281339-.6737426063024382240549898*1j, - -.7250839687106612822281339+.6737426063024382240549898*1j, - -.6506315378609463397807996-.7920349342629491368548074*1j, - -.6506315378609463397807996+.7920349342629491368548074*1j, - -.5564766488918562465935297-.9148198405846724121600860*1j, - -.5564766488918562465935297+.9148198405846724121600860*1j, - -.4345168906815271799687308-1.045382255856986531461592*1j, - -.4345168906815271799687308+1.045382255856986531461592*1j, - -.2640041595834031147954813-1.192762031948052470183960*1j, - -.2640041595834031147954813+1.192762031948052470183960*1j] + -.9025428073192696303995083 - .1105252572789856480992275j, + -.9025428073192696303995083 + .1105252572789856480992275j, + -.8883808106664449854431605 - .2213069215084350419975358j, + -.8883808106664449854431605 + .2213069215084350419975358j, + -.8643915813643204553970169 - .3326258512522187083009453j, + -.8643915813643204553970169 + .3326258512522187083009453j, + -.8299435470674444100273463 - .4448177739407956609694059j, + -.8299435470674444100273463 + .4448177739407956609694059j, + -.7840287980408341576100581 - .5583186348022854707564856j, + -.7840287980408341576100581 + .5583186348022854707564856j, + -.7250839687106612822281339 - .6737426063024382240549898j, + -.7250839687106612822281339 + .6737426063024382240549898j, + -.6506315378609463397807996 - .7920349342629491368548074j, + -.6506315378609463397807996 + .7920349342629491368548074j, + -.5564766488918562465935297 - .9148198405846724121600860j, + -.5564766488918562465935297 + .9148198405846724121600860j, + -.4345168906815271799687308 - 1.045382255856986531461592j, + -.4345168906815271799687308 + 1.045382255856986531461592j, + -.2640041595834031147954813 - 1.192762031948052470183960j, + -.2640041595834031147954813 + 1.192762031948052470183960j] elif N == 22: - p = [-.9058702269930872551848625-52774908289999045189007100.0e-27*1j, - -.9058702269930872551848625+52774908289999045189007100.0e-27*1j, - -.8972983138153530955952835-.1584351912289865608659759*1j, - -.8972983138153530955952835+.1584351912289865608659759*1j, - -.8799661455640176154025352-.2644363039201535049656450*1j, - -.8799661455640176154025352+.2644363039201535049656450*1j, - -.8534754036851687233084587-.3710389319482319823405321*1j, - -.8534754036851687233084587+.3710389319482319823405321*1j, - -.8171682088462720394344996-.4785619492202780899653575*1j, - -.8171682088462720394344996+.4785619492202780899653575*1j, - -.7700332930556816872932937-.5874255426351153211965601*1j, - -.7700332930556816872932937+.5874255426351153211965601*1j, - -.7105305456418785989070935-.6982266265924524000098548*1j, - -.7105305456418785989070935+.6982266265924524000098548*1j, - -.6362427683267827226840153-.8118875040246347267248508*1j, - -.6362427683267827226840153+.8118875040246347267248508*1j, - -.5430983056306302779658129-.9299947824439872998916657*1j, - -.5430983056306302779658129+.9299947824439872998916657*1j, - -.4232528745642628461715044-1.055755605227545931204656*1j, - -.4232528745642628461715044+1.055755605227545931204656*1j, - -.2566376987939318038016012-1.197982433555213008346532*1j, - -.2566376987939318038016012+1.197982433555213008346532*1j] + p = [-.9058702269930872551848625 - 52774908289999045189007100.0e-27j, + -.9058702269930872551848625 + 52774908289999045189007100.0e-27j, + -.8972983138153530955952835 - .1584351912289865608659759j, + -.8972983138153530955952835 + .1584351912289865608659759j, + -.8799661455640176154025352 - .2644363039201535049656450j, + -.8799661455640176154025352 + .2644363039201535049656450j, + -.8534754036851687233084587 - .3710389319482319823405321j, + -.8534754036851687233084587 + .3710389319482319823405321j, + -.8171682088462720394344996 - .4785619492202780899653575j, + -.8171682088462720394344996 + .4785619492202780899653575j, + -.7700332930556816872932937 - .5874255426351153211965601j, + -.7700332930556816872932937 + .5874255426351153211965601j, + -.7105305456418785989070935 - .6982266265924524000098548j, + -.7105305456418785989070935 + .6982266265924524000098548j, + -.6362427683267827226840153 - .8118875040246347267248508j, + -.6362427683267827226840153 + .8118875040246347267248508j, + -.5430983056306302779658129 - .9299947824439872998916657j, + -.5430983056306302779658129 + .9299947824439872998916657j, + -.4232528745642628461715044 - 1.055755605227545931204656j, + -.4232528745642628461715044 + 1.055755605227545931204656j, + -.2566376987939318038016012 - 1.197982433555213008346532j, + -.2566376987939318038016012 + 1.197982433555213008346532j] elif N == 23: p = [-.9066732476324988168207439, - -.9027564979912504609412993-.1010534335314045013252480*1j, - -.9027564979912504609412993+.1010534335314045013252480*1j, - -.8909283242471251458653994-.2023024699381223418195228*1j, - -.8909283242471251458653994+.2023024699381223418195228*1j, - -.8709469395587416239596874-.3039581993950041588888925*1j, - -.8709469395587416239596874+.3039581993950041588888925*1j, - -.8423805948021127057054288-.4062657948237602726779246*1j, - -.8423805948021127057054288+.4062657948237602726779246*1j, - -.8045561642053176205623187-.5095305912227258268309528*1j, - -.8045561642053176205623187+.5095305912227258268309528*1j, - -.7564660146829880581478138-.6141594859476032127216463*1j, - -.7564660146829880581478138+.6141594859476032127216463*1j, - -.6965966033912705387505040-.7207341374753046970247055*1j, - -.6965966033912705387505040+.7207341374753046970247055*1j, - -.6225903228771341778273152-.8301558302812980678845563*1j, - -.6225903228771341778273152+.8301558302812980678845563*1j, - -.5304922463810191698502226-.9439760364018300083750242*1j, - -.5304922463810191698502226+.9439760364018300083750242*1j, - -.4126986617510148836149955-1.065328794475513585531053*1j, - -.4126986617510148836149955+1.065328794475513585531053*1j, - -.2497697202208956030229911-1.202813187870697831365338*1j, - -.2497697202208956030229911+1.202813187870697831365338*1j] + -.9027564979912504609412993 - .1010534335314045013252480j, + -.9027564979912504609412993 + .1010534335314045013252480j, + -.8909283242471251458653994 - .2023024699381223418195228j, + -.8909283242471251458653994 + .2023024699381223418195228j, + -.8709469395587416239596874 - .3039581993950041588888925j, + -.8709469395587416239596874 + .3039581993950041588888925j, + -.8423805948021127057054288 - .4062657948237602726779246j, + -.8423805948021127057054288 + .4062657948237602726779246j, + -.8045561642053176205623187 - .5095305912227258268309528j, + -.8045561642053176205623187 + .5095305912227258268309528j, + -.7564660146829880581478138 - .6141594859476032127216463j, + -.7564660146829880581478138 + .6141594859476032127216463j, + -.6965966033912705387505040 - .7207341374753046970247055j, + -.6965966033912705387505040 + .7207341374753046970247055j, + -.6225903228771341778273152 - .8301558302812980678845563j, + -.6225903228771341778273152 + .8301558302812980678845563j, + -.5304922463810191698502226 - .9439760364018300083750242j, + -.5304922463810191698502226 + .9439760364018300083750242j, + -.4126986617510148836149955 - 1.065328794475513585531053j, + -.4126986617510148836149955 + 1.065328794475513585531053j, + -.2497697202208956030229911 - 1.202813187870697831365338j, + -.2497697202208956030229911 + 1.202813187870697831365338j] elif N == 24: - p = [-.9055312363372773709269407-48440066540478700874836350.0e-27*1j, - -.9055312363372773709269407+48440066540478700874836350.0e-27*1j, - -.8983105104397872954053307-.1454056133873610120105857*1j, - -.8983105104397872954053307+.1454056133873610120105857*1j, - -.8837358034555706623131950-.2426335234401383076544239*1j, - -.8837358034555706623131950+.2426335234401383076544239*1j, - -.8615278304016353651120610-.3403202112618624773397257*1j, - -.8615278304016353651120610+.3403202112618624773397257*1j, - -.8312326466813240652679563-.4386985933597305434577492*1j, - -.8312326466813240652679563+.4386985933597305434577492*1j, - -.7921695462343492518845446-.5380628490968016700338001*1j, - -.7921695462343492518845446+.5380628490968016700338001*1j, - -.7433392285088529449175873-.6388084216222567930378296*1j, - -.7433392285088529449175873+.6388084216222567930378296*1j, - -.6832565803536521302816011-.7415032695091650806797753*1j, - -.6832565803536521302816011+.7415032695091650806797753*1j, - -.6096221567378335562589532-.8470292433077202380020454*1j, - -.6096221567378335562589532+.8470292433077202380020454*1j, - -.5185914574820317343536707-.9569048385259054576937721*1j, - -.5185914574820317343536707+.9569048385259054576937721*1j, - -.4027853855197518014786978-1.074195196518674765143729*1j, - -.4027853855197518014786978+1.074195196518674765143729*1j, - -.2433481337524869675825448-1.207298683731972524975429*1j, - -.2433481337524869675825448+1.207298683731972524975429*1j] + p = [-.9055312363372773709269407 - 48440066540478700874836350.0e-27j, + -.9055312363372773709269407 + 48440066540478700874836350.0e-27j, + -.8983105104397872954053307 - .1454056133873610120105857j, + -.8983105104397872954053307 + .1454056133873610120105857j, + -.8837358034555706623131950 - .2426335234401383076544239j, + -.8837358034555706623131950 + .2426335234401383076544239j, + -.8615278304016353651120610 - .3403202112618624773397257j, + -.8615278304016353651120610 + .3403202112618624773397257j, + -.8312326466813240652679563 - .4386985933597305434577492j, + -.8312326466813240652679563 + .4386985933597305434577492j, + -.7921695462343492518845446 - .5380628490968016700338001j, + -.7921695462343492518845446 + .5380628490968016700338001j, + -.7433392285088529449175873 - .6388084216222567930378296j, + -.7433392285088529449175873 + .6388084216222567930378296j, + -.6832565803536521302816011 - .7415032695091650806797753j, + -.6832565803536521302816011 + .7415032695091650806797753j, + -.6096221567378335562589532 - .8470292433077202380020454j, + -.6096221567378335562589532 + .8470292433077202380020454j, + -.5185914574820317343536707 - .9569048385259054576937721j, + -.5185914574820317343536707 + .9569048385259054576937721j, + -.4027853855197518014786978 - 1.074195196518674765143729j, + -.4027853855197518014786978 + 1.074195196518674765143729j, + -.2433481337524869675825448 - 1.207298683731972524975429j, + -.2433481337524869675825448 + 1.207298683731972524975429j] elif N == 25: p = [-.9062073871811708652496104, - -.9028833390228020537142561-93077131185102967450643820.0e-27*1j, - -.9028833390228020537142561+93077131185102967450643820.0e-27*1j, - -.8928551459883548836774529-.1863068969804300712287138*1j, - -.8928551459883548836774529+.1863068969804300712287138*1j, - -.8759497989677857803656239-.2798521321771408719327250*1j, - -.8759497989677857803656239+.2798521321771408719327250*1j, - -.8518616886554019782346493-.3738977875907595009446142*1j, - -.8518616886554019782346493+.3738977875907595009446142*1j, - -.8201226043936880253962552-.4686668574656966589020580*1j, - -.8201226043936880253962552+.4686668574656966589020580*1j, - -.7800496278186497225905443-.5644441210349710332887354*1j, - -.7800496278186497225905443+.5644441210349710332887354*1j, - -.7306549271849967721596735-.6616149647357748681460822*1j, - -.7306549271849967721596735+.6616149647357748681460822*1j, - -.6704827128029559528610523-.7607348858167839877987008*1j, - -.6704827128029559528610523+.7607348858167839877987008*1j, - -.5972898661335557242320528-.8626676330388028512598538*1j, - -.5972898661335557242320528+.8626676330388028512598538*1j, - -.5073362861078468845461362-.9689006305344868494672405*1j, - -.5073362861078468845461362+.9689006305344868494672405*1j, - -.3934529878191079606023847-1.082433927173831581956863*1j, - -.3934529878191079606023847+1.082433927173831581956863*1j, - -.2373280669322028974199184-1.211476658382565356579418*1j, - -.2373280669322028974199184+1.211476658382565356579418*1j] + -.9028833390228020537142561 - 93077131185102967450643820.0e-27j, + -.9028833390228020537142561 + 93077131185102967450643820.0e-27j, + -.8928551459883548836774529 - .1863068969804300712287138j, + -.8928551459883548836774529 + .1863068969804300712287138j, + -.8759497989677857803656239 - .2798521321771408719327250j, + -.8759497989677857803656239 + .2798521321771408719327250j, + -.8518616886554019782346493 - .3738977875907595009446142j, + -.8518616886554019782346493 + .3738977875907595009446142j, + -.8201226043936880253962552 - .4686668574656966589020580j, + -.8201226043936880253962552 + .4686668574656966589020580j, + -.7800496278186497225905443 - .5644441210349710332887354j, + -.7800496278186497225905443 + .5644441210349710332887354j, + -.7306549271849967721596735 - .6616149647357748681460822j, + -.7306549271849967721596735 + .6616149647357748681460822j, + -.6704827128029559528610523 - .7607348858167839877987008j, + -.6704827128029559528610523 + .7607348858167839877987008j, + -.5972898661335557242320528 - .8626676330388028512598538j, + -.5972898661335557242320528 + .8626676330388028512598538j, + -.5073362861078468845461362 - .9689006305344868494672405j, + -.5073362861078468845461362 + .9689006305344868494672405j, + -.3934529878191079606023847 - 1.082433927173831581956863j, + -.3934529878191079606023847 + 1.082433927173831581956863j, + -.2373280669322028974199184 - 1.211476658382565356579418j, + -.2373280669322028974199184 + 1.211476658382565356579418j] else: raise ValueError("Bessel Filter not supported for order %d" % N) return z, p, k -filter_dict = {'butter': [buttap,buttord], - 'butterworth' : [buttap,buttord], - 'cauer' : [ellipap,ellipord], - 'elliptic' : [ellipap,ellipord], - 'ellip' : [ellipap,ellipord], - 'bessel' : [besselap], - 'cheby1' : [cheb1ap, cheb1ord], - 'chebyshev1' : [cheb1ap, cheb1ord], - 'chebyshevi' : [cheb1ap, cheb1ord], - 'cheby2' : [cheb2ap, cheb2ord], - 'chebyshev2' : [cheb2ap, cheb2ord], - 'chebyshevii' : [cheb2ap, cheb2ord] +filter_dict = {'butter': [buttap, buttord], + 'butterworth': [buttap, buttord], + 'cauer': [ellipap, ellipord], + 'elliptic': [ellipap, ellipord], + 'ellip': [ellipap, ellipord], + 'bessel': [besselap], + 'cheby1': [cheb1ap, cheb1ord], + 'chebyshev1': [cheb1ap, cheb1ord], + 'chebyshevi': [cheb1ap, cheb1ord], + 'cheby2': [cheb2ap, cheb2ord], + 'chebyshev2': [cheb2ap, cheb2ord], + 'chebyshevii': [cheb2ap, cheb2ord] } -band_dict = {'band':'bandpass', - 'bandpass':'bandpass', - 'pass' : 'bandpass', - 'bp':'bandpass', - 'bs':'bandstop', - 'bandstop':'bandstop', - 'bands' : 'bandstop', - 'stop' : 'bandstop', - 'l' : 'lowpass', +band_dict = {'band': 'bandpass', + 'bandpass': 'bandpass', + 'pass': 'bandpass', + 'bp': 'bandpass', + 'bs': 'bandstop', + 'bandstop': 'bandstop', + 'bands': 'bandstop', + 'stop': 'bandstop', + 'l': 'lowpass', 'low': 'lowpass', - 'lowpass' : 'lowpass', - 'high' : 'highpass', - 'highpass' : 'highpass', - 'h' : 'highpass' + 'lowpass': 'lowpass', + 'high': 'highpass', + 'highpass': 'highpass', + 'h': 'highpass' } warnings.simplefilter("always", BadCoefficients) diff --git a/scipy/signal/fir_filter_design.py b/scipy/signal/fir_filter_design.py index db5467148b2c..405f3f21f32f 100644 --- a/scipy/signal/fir_filter_design.py +++ b/scipy/signal/fir_filter_design.py @@ -46,7 +46,7 @@ def kaiser_beta(a): if a > 50: beta = 0.1102 * (a - 8.7) elif a > 21: - beta = 0.5842 * (a - 21)**0.4 + 0.07886 * (a - 21) + beta = 0.5842 * (a - 21) ** 0.4 + 0.07886 * (a - 21) else: beta = 0.0 return beta @@ -65,8 +65,8 @@ def kaiser_atten(numtaps, width): N : int The number of taps in the FIR filter. width : float - The desired width of the transition region between passband and stopband - (or, in general, at any discontinuity) for the filter. + The desired width of the transition region between passband and + stopband (or, in general, at any discontinuity) for the filter. Returns ------- @@ -138,9 +138,9 @@ def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, """ FIR filter design using the window method. - This function computes the coefficients of a finite impulse response filter. - The filter will have linear phase; it will be Type I if `numtaps` is odd and - Type II if `numtaps` is even. + This function computes the coefficients of a finite impulse response + filter. The filter will have linear phase; it will be Type I if + `numtaps` is odd and Type II if `numtaps` is even. Type II filters always have zero response at the Nyquist rate, so a ValueError exception is raised if firwin is called with `numtaps` even and @@ -245,18 +245,21 @@ def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, # Check for invalid input. if cutoff.ndim > 1: - raise ValueError("The cutoff argument must be at most one-dimensional.") + raise ValueError("The cutoff argument must be at most " + "one-dimensional.") if cutoff.size == 0: raise ValueError("At least one cutoff frequency must be given.") if cutoff.min() <= 0 or cutoff.max() >= 1: - raise ValueError("Invalid cutoff frequency: frequencies must be greater than 0 and less than nyq.") + raise ValueError("Invalid cutoff frequency: frequencies must be " + "greater than 0 and less than nyq.") if np.any(np.diff(cutoff) <= 0): - raise ValueError("Invalid cutoff frequencies: the frequencies must be strictly increasing.") + raise ValueError("Invalid cutoff frequencies: the frequencies " + "must be strictly increasing.") if width is not None: # A width was given. Find the beta parameter of the Kaiser window # and set `window`. This overrides the value of `window` passed in. - atten = kaiser_atten(numtaps, float(width)/nyq) + atten = kaiser_atten(numtaps, float(width) / nyq) beta = kaiser_beta(atten) window = ('kaiser', beta) @@ -265,15 +268,16 @@ def firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, raise ValueError("A filter with an even number of coefficients must " "have zero response at the Nyquist rate.") - # Insert 0 and/or 1 at the ends of cutoff so that the length of cutoff is even, - # and each pair in cutoff corresponds to passband. - cutoff = np.hstack(([0.0]*pass_zero, cutoff, [1.0]*pass_nyquist)) + # Insert 0 and/or 1 at the ends of cutoff so that the length of cutoff + # is even, and each pair in cutoff corresponds to passband. + cutoff = np.hstack(([0.0] * pass_zero, cutoff, [1.0] * pass_nyquist)) - # `bands` is a 2D array; each row gives the left and right edges of a passband. - bands = cutoff.reshape(-1,2) + # `bands` is a 2D array; each row gives the left and right edges of + # a passband. + bands = cutoff.reshape(-1, 2) # Build up the coefficients. - alpha = 0.5 * (numtaps-1) + alpha = 0.5 * (numtaps - 1) m = np.arange(0, numtaps) - alpha h = 0 for left, right in bands: @@ -397,8 +401,9 @@ def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=1.0): raise ValueError('freq and gain must be of same length.') if nfreqs is not None and numtaps >= nfreqs: - raise ValueError('ntaps must be less than nfreqs, but firwin2 was ' - 'called with ntaps=%d and nfreqs=%s' % (numtaps, nfreqs)) + raise ValueError(('ntaps must be less than nfreqs, but firwin2 was ' + 'called with ntaps=%d and nfreqs=%s') % + (numtaps, nfreqs)) if freq[0] != 0 or freq[-1] != nyq: raise ValueError('freq must start with 0 and end with `nyq`.') @@ -414,14 +419,14 @@ def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=1.0): "have zero gain at the Nyquist rate.") if nfreqs is None: - nfreqs = 1 + 2 ** int(ceil(log(numtaps,2))) + nfreqs = 1 + 2 ** int(ceil(log(numtaps, 2))) # Tweak any repeated values in freq so that interp works. eps = np.finfo(float).eps for k in range(len(freq)): - if k < len(freq)-1 and freq[k] == freq[k+1]: + if k < len(freq) - 1 and freq[k] == freq[k + 1]: freq[k] = freq[k] - eps - freq[k+1] = freq[k+1] + eps + freq[k + 1] = freq[k + 1] + eps # Linearly interpolate the desired response on a uniform mesh `x`. x = np.linspace(0.0, nyq, nfreqs) @@ -429,7 +434,7 @@ def firwin2(numtaps, freq, gain, nfreqs=None, window='hamming', nyq=1.0): # Adjust the phases of the coefficients so that the first `ntaps` of the # inverse FFT are the desired filter coefficients. - shift = np.exp(-(numtaps-1)/2. * 1.j * np.pi * x / nyq) + shift = np.exp(-(numtaps - 1) / 2. * 1.j * np.pi * x / nyq) fx2 = fx * shift # Use irfft to compute the inverse FFT. @@ -534,9 +539,10 @@ def remez(numtaps, bands, desired, weight=None, Hz=1, type='bandpass', """ # Convert type try: - tnum = {'bandpass':1, 'differentiator':2, 'hilbert':3}[type] + tnum = {'bandpass': 1, 'differentiator': 2, 'hilbert': 3}[type] except KeyError: - raise ValueError("Type must be 'bandpass', 'differentiator', or 'hilbert'") + raise ValueError("Type must be 'bandpass', 'differentiator', " + "or 'hilbert'") # Convert weight if weight is None: diff --git a/scipy/signal/info.py b/scipy/signal/info.py index 740a83cb764b..344ee00d49cd 100644 --- a/scipy/signal/info.py +++ b/scipy/signal/info.py @@ -69,9 +69,11 @@ Filter design: bilinear: - Return a digital filter from an analog filter using the bilinear transform. + Return a digital filter from an analog filter using the bilinear + transform. firwin: - Windowed FIR filter design, with frequency response defined as pass and stop bands. + Windowed FIR filter design, with frequency response defined as + pass and stop bands. firwin2: Windowed FIR filter design, with arbitrary frequency response. freqs: @@ -85,10 +87,12 @@ invres: Inverse partial fraction expansion. kaiser_beta: - Compute the Kaiser parameter beta, given the desired FIR filter attenuation. + Compute the Kaiser parameter beta, given the desired FIR filter + attenuation. kaiser_atten: - Compute the attenuation of a Kaiser FIR filter, given the number of taps - and the transition width at discontinuities in the frequency response. + Compute the attenuation of a Kaiser FIR filter, given the number of + taps and the transition width at discontinuities in the frequency + response. kaiserord: Design a Kaiser window to limit ripple and width of transition region. remez: diff --git a/scipy/signal/ltisys.py b/scipy/signal/ltisys.py index bc42c6e557be..eb2cfce61112 100644 --- a/scipy/signal/ltisys.py +++ b/scipy/signal/ltisys.py @@ -52,35 +52,37 @@ def tf2ss(num, den): num = asarray([num], num.dtype) M = num.shape[1] K = len(den) - if (M > K): + if M > K: raise ValueError("Improper transfer function.") - if (M == 0 or K == 0): # Null system - return array([],float), array([], float), array([], float), \ + if M == 0 or K == 0: # Null system + return array([], float), array([], float), array([], float), \ array([], float) # pad numerator to have same number of columns has denominator - num = r_['-1',zeros((num.shape[0],K-M), num.dtype), num] + num = r_['-1', zeros((num.shape[0], K - M), num.dtype), num] if num.shape[-1] > 0: - D = num[:,0] + D = num[:, 0] else: - D = array([],float) + D = array([], float) if K == 1: return array([], float), array([], float), array([], float), D frow = -array([den[1:]]) - A = r_[frow, eye(K-2, K-1)] - B = eye(K-1, 1) - C = num[:,1:] - num[:,0] * den[1:] + A = r_[frow, eye(K - 2, K - 1)] + B = eye(K - 1, 1) + C = num[:, 1:] - num[:, 0] * den[1:] return A, B, C, D + def _none_to_empty(arg): if arg is None: return [] else: return arg + def abcd_normalize(A=None, B=None, C=None, D=None): """Check state-space matrices and ensure they are rank-2. @@ -123,6 +125,7 @@ def abcd_normalize(A=None, B=None, C=None, D=None): return A, B, C, D + def ss2tf(A, B, C, D, input=0): """State-space to transfer function. @@ -152,31 +155,32 @@ def ss2tf(A, B, C, D, input=0): # make MOSI from possibly MOMI system. if B.shape[-1] != 0: - B = B[:,input] - B.shape = (B.shape[0],1) + B = B[:, input] + B.shape = (B.shape[0], 1) if D.shape[-1] != 0: - D = D[:,input] + D = D[:, input] try: den = poly(A) except ValueError: den = 1 - if (product(B.shape,axis=0) == 0) and (product(C.shape,axis=0) == 0): + if (product(B.shape, axis=0) == 0) and (product(C.shape, axis=0) == 0): num = numpy.ravel(D) - if (product(D.shape,axis=0) == 0) and (product(A.shape,axis=0) == 0): + if (product(D.shape, axis=0) == 0) and (product(A.shape, axis=0) == 0): den = [] return num, den num_states = A.shape[0] - type_test = A[:,0] + B[:,0] + C[0,:] + D - num = numpy.zeros((nout, num_states+1), type_test.dtype) + type_test = A[:, 0] + B[:, 0] + C[0, :] + D + num = numpy.zeros((nout, num_states + 1), type_test.dtype) for k in range(nout): - Ck = atleast_2d(C[k,:]) - num[k] = poly(A - dot(B,Ck)) + (D[k]-1)*den + Ck = atleast_2d(C[k, :]) + num[k] = poly(A - dot(B, Ck)) + (D[k] - 1) * den return num, den + def zpk2ss(z, p, k): """Zero-pole-gain representation to state-space representation @@ -193,7 +197,8 @@ def zpk2ss(z, p, k): State-space matrices. """ - return tf2ss(*zpk2tf(z,p,k)) + return tf2ss(*zpk2tf(z, p, k)) + def ss2zpk(A, B, C, D, input=0): """State-space representation to zero-pole-gain representation. @@ -213,12 +218,13 @@ def ss2zpk(A, B, C, D, input=0): System gain. """ - return tf2zpk(*ss2tf(A,B,C,D,input=input)) + return tf2zpk(*ss2tf(A, B, C, D, input=input)) + class lti(object): """Linear Time Invariant class which simplifies representation. """ - def __init__(self,*args,**kwords): + def __init__(self, *args, **kwords): """Initialize the LTI system using either: (numerator, denominator) (zeros, poles, gain) @@ -262,7 +268,7 @@ def __init__(self,*args,**kwords): raise ValueError("Needs 2, 3, or 4 arguments.") def __setattr__(self, attr, val): - if attr in ['num','den']: + if attr in ['num', 'den']: self.__dict__[attr] = val self.__dict__['zeros'], self.__dict__['poles'], \ self.__dict__['gain'] = \ @@ -358,7 +364,7 @@ def lsim2(system, U=None, T=None, X0=None, **kwargs): sys = lti(*system) if X0 is None: - X0 = zeros(sys.B.shape[0],sys.A.dtype) + X0 = zeros(sys.B.shape[0], sys.A.dtype) if T is None: # XXX T should really be a required argument, but U was @@ -376,7 +382,7 @@ def lsim2(system, U=None, T=None, X0=None, **kwargs): if U is not None: U = atleast_1d(U) if len(U.shape) == 1: - U = U.reshape(-1,1) + U = U.reshape(-1, 1) sU = U.shape if sU[0] != len(T): raise ValueError("U must have the same number of rows " @@ -393,15 +399,15 @@ def lsim2(system, U=None, T=None, X0=None, **kwargs): def fprime(x, t, sys, ufunc): """The vector field of the linear system.""" - return dot(sys.A,x) + squeeze(dot(sys.B,nan_to_num(ufunc([t])))) + return dot(sys.A, x) + squeeze(dot(sys.B, nan_to_num(ufunc([t])))) xout = integrate.odeint(fprime, X0, T, args=(sys, ufunc), **kwargs) - yout = dot(sys.C,transpose(xout)) + dot(sys.D,transpose(U)) + yout = dot(sys.C, transpose(xout)) + dot(sys.D, transpose(U)) else: def fprime(x, t, sys): """The vector field of the linear system.""" - return dot(sys.A,x) + return dot(sys.A, x) xout = integrate.odeint(fprime, X0, T, args=(sys,), **kwargs) - yout = dot(sys.C,transpose(xout)) + yout = dot(sys.C, transpose(xout)) return T, squeeze(transpose(yout)), xout @@ -459,7 +465,7 @@ def lsim(system, U, T, X0=None, interp=1): U = atleast_1d(U) T = atleast_1d(T) if len(U.shape) == 1: - U = U.reshape((U.shape[0],1)) + U = U.reshape((U.shape[0], 1)) sU = U.shape if len(T.shape) != 1: raise ValueError("T must be a rank-1 array.") @@ -472,38 +478,40 @@ def lsim(system, U, T, X0=None, interp=1): if X0 is None: X0 = zeros(sys.B.shape[0], sys.A.dtype) - xout = zeros((len(T),sys.B.shape[0]), sys.A.dtype) + xout = zeros((len(T), sys.B.shape[0]), sys.A.dtype) xout[0] = X0 A = sys.A AT, BT = transpose(sys.A), transpose(sys.B) - dt = T[1]-T[0] + dt = T[1] - T[0] lam, v = linalg.eig(A) vt = transpose(v) vti = linalg.inv(vt) - GT = dot(dot(vti,diag(numpy.exp(dt*lam))),vt).astype(xout.dtype) + GT = dot(dot(vti, diag(numpy.exp(dt * lam))), vt).astype(xout.dtype) ATm1 = linalg.inv(AT) - ATm2 = dot(ATm1,ATm1) - I = eye(A.shape[0],dtype=A.dtype) - GTmI = GT-I - F1T = dot(dot(BT,GTmI),ATm1) + ATm2 = dot(ATm1, ATm1) + I = eye(A.shape[0], dtype=A.dtype) + GTmI = GT - I + F1T = dot(dot(BT, GTmI), ATm1) if interp: - F2T = dot(BT,dot(GTmI,ATm2)/dt - ATm1) + F2T = dot(BT, dot(GTmI, ATm2) / dt - ATm1) - for k in xrange(1,len(T)): - dt1 = T[k] - T[k-1] + for k in xrange(1, len(T)): + dt1 = T[k] - T[k - 1] if dt1 != dt: dt = dt1 - GT = dot(dot(vti,diag(numpy.exp(dt*lam))),vt).astype(xout.dtype) - GTmI = GT-I - F1T = dot(dot(BT,GTmI),ATm1) + GT = dot(dot(vti, diag(numpy.exp(dt * lam))), + vt).astype(xout.dtype) + GTmI = GT - I + F1T = dot(dot(BT, GTmI), ATm1) if interp: - F2T = dot(BT,dot(GTmI,ATm2)/dt - ATm1) + F2T = dot(BT, dot(GTmI, ATm2) / dt - ATm1) - xout[k] = dot(xout[k-1],GT) + dot(U[k-1],F1T) + xout[k] = dot(xout[k - 1], GT) + dot(U[k - 1], F1T) if interp: - xout[k] = xout[k] + dot((U[k]-U[k-1]),F2T) + xout[k] = xout[k] + dot((U[k] - U[k - 1]), F2T) - yout = squeeze(dot(U,transpose(sys.D))) + squeeze(dot(xout,transpose(sys.C))) + yout = (squeeze(dot(U, transpose(sys.D))) + + squeeze(dot(xout, transpose(sys.C)))) return T, squeeze(yout), squeeze(xout) @@ -534,7 +542,7 @@ def _default_response_times(A, n): if r == 0.0: r = 1.0 tc = 1.0 / r - t = linspace(0.0, 7*tc, n) + t = linspace(0.0, 7 * tc, n) return t @@ -575,13 +583,13 @@ def impulse(system, X0=None, T=None, N=None): if T is None: T = _default_response_times(sys.A, N) h = zeros(T.shape, sys.A.dtype) - s,v = linalg.eig(sys.A) + s, v = linalg.eig(sys.A) vi = linalg.inv(v) C = sys.C for k in range(len(h)): - es = diag(numpy.exp(s*T[k])) - eA = (dot(dot(v,es),vi)).astype(h.dtype) - h[k] = squeeze(dot(dot(C,eA),B)) + es = diag(numpy.exp(s * T[k])) + eA = (dot(dot(v, es), vi)).astype(h.dtype) + h[k] = squeeze(dot(dot(C, eA), B)) return T, h @@ -707,6 +715,7 @@ def step(system, X0=None, T=None, N=None): vals = lsim(sys, U, T, X0=X0) return vals[0], vals[1] + def step2(system, X0=None, T=None, N=None, **kwargs): """Step response of continuous-time system. diff --git a/scipy/signal/setup.py b/scipy/signal/setup.py index 0064c545e988..34f84e30cd62 100755 --- a/scipy/signal/setup.py +++ b/scipy/signal/setup.py @@ -1,6 +1,7 @@ #!/usr/bin/env python -def configuration(parent_package='',top_path=None): + +def configuration(parent_package='', top_path=None): from numpy.distutils.misc_util import Configuration config = Configuration('signal', parent_package, top_path) @@ -8,22 +9,23 @@ def configuration(parent_package='',top_path=None): config.add_data_dir('tests') config.add_extension('sigtools', - sources=['sigtoolsmodule.c', - 'firfilter.c','medianfilter.c', 'lfilter.c.src', + sources=['sigtoolsmodule.c', 'firfilter.c', + 'medianfilter.c', 'lfilter.c.src', 'correlate_nd.c.src'], - depends = ['sigtools.h'], + depends=['sigtools.h'], include_dirs=['.'] ) config.add_extension('spectral', sources=['spectral.c']) config.add_extension('spline', - sources = ['splinemodule.c','S_bspline_util.c','D_bspline_util.c', - 'C_bspline_util.c','Z_bspline_util.c','bspline_util.c'], + sources=['splinemodule.c', 'S_bspline_util.c', 'D_bspline_util.c', + 'C_bspline_util.c', 'Z_bspline_util.c', 'bspline_util.c'], ) return config + if __name__ == '__main__': from numpy.distutils.core import setup setup(**configuration(top_path='').todict()) diff --git a/scipy/signal/setupscons.py b/scipy/signal/setupscons.py index 21488955b7df..0a3009c2f526 100755 --- a/scipy/signal/setupscons.py +++ b/scipy/signal/setupscons.py @@ -1,6 +1,7 @@ #!/usr/bin/env python -def configuration(parent_package='',top_path=None): + +def configuration(parent_package='', top_path=None): from numpy.distutils.misc_util import Configuration config = Configuration('signal', parent_package, top_path) @@ -10,6 +11,7 @@ def configuration(parent_package='',top_path=None): return config + if __name__ == '__main__': from numpy.distutils.core import setup setup(**configuration(top_path='').todict()) diff --git a/scipy/signal/signaltools.py b/scipy/signal/signaltools.py index f2e342a69267..11dbc05d54ee 100644 --- a/scipy/signal/signaltools.py +++ b/scipy/signal/signaltools.py @@ -8,7 +8,7 @@ ifftn, fftfreq from numpy import polyadd, polymul, polydiv, polysub, roots, \ poly, polyval, polyder, cast, asarray, isscalar, atleast_1d, \ - ones, real, real_if_close, zeros, array, arange, where, rank, \ + ones, real_if_close, zeros, array, arange, where, rank, \ newaxis, product, ravel, sum, r_, iscomplexobj, take, \ argsort, allclose, expand_dims, unique, prod, sort, reshape, \ transpose, dot, mean, flipud, ndarray, atleast_2d @@ -735,7 +735,7 @@ def hilbert2(x, N=None): if len(x.shape) > 2: raise ValueError("x must be rank 2.") if iscomplexobj(x): - raise ValueError("x must be real.") + raise ValueError("x must be real.") if N is None: N = x.shape elif isinstance(N, int): diff --git a/scipy/signal/waveforms.py b/scipy/signal/waveforms.py index bf59df2d30f2..30597a67b0dd 100644 --- a/scipy/signal/waveforms.py +++ b/scipy/signal/waveforms.py @@ -124,7 +124,8 @@ def square(t, duty=0.5): return y -def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False, retenv=False): +def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False, + retenv=False): """ Return a gaussian modulated sinusoid: exp(-a t^2) exp(1j*2*pi*fc*t). @@ -310,8 +311,8 @@ def _chirp_phase(t, f0, t1, f1, method='linear', vertex_zero=True): elif method in ['logarithmic', 'log', 'lo']: if f0 * f1 <= 0.0: - raise ValueError("For a geometric chirp, f0 and f1 must be nonzero " \ - "and have the same sign.") + raise ValueError("For a geometric chirp, f0 and f1 must be " + "nonzero and have the same sign.") if f0 == f1: phase = 2 * pi * f0 * t else: @@ -326,8 +327,8 @@ def _chirp_phase(t, f0, t1, f1, method='linear', vertex_zero=True): phase = 2 * pi * (f0 * c / df) * log((df * t + c) / c) else: - raise ValueError("method must be 'linear', 'quadratic', 'logarithmic', " - "or 'hyperbolic', but a value of %r was given." % method) + raise ValueError("method must be 'linear', 'quadratic', 'logarithmic'," + " or 'hyperbolic', but a value of %r was given." % method) return phase diff --git a/scipy/signal/wavelets.py b/scipy/signal/wavelets.py index dd220350547a..3db241215453 100644 --- a/scipy/signal/wavelets.py +++ b/scipy/signal/wavelets.py @@ -28,32 +28,34 @@ def daub(p): elif p == 2: f = sqrt(2) / 8 c = sqrt(3) - return f * np.array([1+c, 3+c, 3-c, 1-c]) + return f * np.array([1 + c, 3 + c, 3 - c, 1 - c]) elif p == 3: tmp = 12 * sqrt(10) - z1 = 1.5 + sqrt(15+tmp)/6 - 1j*(sqrt(15)+sqrt(tmp-15))/6 + z1 = 1.5 + sqrt(15 + tmp) / 6 - 1j * (sqrt(15) + sqrt(tmp - 15)) / 6 z1c = np.conj(z1) - f = sqrt(2)/8 - d0 = np.real((1-z1)*(1-z1c)) - a0 = np.real(z1*z1c) - a1 = 2*np.real(z1) - return f/d0*np.array([a0, 3*a0-a1, 3*a0-3*a1+1, a0-3*a1+3, 3-a1, 1]) + f = sqrt(2) / 8 + d0 = np.real((1 - z1) * (1 - z1c)) + a0 = np.real(z1 * z1c) + a1 = 2 * np.real(z1) + return f / d0 * np.array([a0, 3 * a0 - a1, 3 * a0 - 3 * a1 + 1, + a0 - 3 * a1 + 3, 3 - a1, 1]) elif p < 35: # construct polynomial and factor it if p < 35: - P = [comb(p-1+k, k, exact=1) for k in range(p)][::-1] + P = [comb(p - 1 + k, k, exact=1) for k in range(p)][::-1] yj = np.roots(P) else: # try different polynomial --- needs work - P = [comb(p-1+k, k, exact=1)/4.0**k for k in range(p)][::-1] + P = [comb(p - 1 + k, k, exact=1) / 4.0 ** k + for k in range(p)][::-1] yj = np.roots(P) / 4 # for each root, compute two z roots, select the one with |z|>1 # Build up final polynomial - c = np.poly1d([1, 1])**p + c = np.poly1d([1, 1]) ** p q = np.poly1d([1]) - for k in range(p-1): + for k in range(p - 1): yval = yj[k] - part = 2*sqrt(yval*(yval-1)) - const = 1 - 2*yval + part = 2 * sqrt(yval * (yval - 1)) + const = 1 - 2 * yval z1 = const + part if (abs(z1)) < 1: z1 = const - part @@ -72,8 +74,8 @@ def qmf(hk): """Return high-pass qmf filter from low-pass """ N = len(hk) - 1 - asgn = [{0:1, 1:-1}[k%2] for k in range(N+1)] - return hk[::-1]*np.array(asgn) + asgn = [{0: 1, 1: -1}[k % 2] for k in range(N + 1)] + return hk[::-1] * np.array(asgn) def wavedec(amn, hk): @@ -122,9 +124,9 @@ def cascade(hk, J=7): """ - N = len(hk)-1 + N = len(hk) - 1 - if (J > 30 - np.log2(N+1)): + if (J > 30 - np.log2(N + 1)): raise ValueError("Too many levels.") if (J < 1): raise ValueError("Too few levels.") @@ -137,8 +139,8 @@ def cascade(hk, J=7): gk = qmf(hk) tgk = np.r_[gk, 0] - indx1 = np.clip(2*nn-kk, -1, N+1) - indx2 = np.clip(2*nn-kk+1, -1, N+1) + indx1 = np.clip(2 * nn - kk, -1, N + 1) + indx2 = np.clip(2 * nn - kk + 1, -1, N + 1) m = np.zeros((2, 2, N, N), 'd') m[0, 0] = np.take(thk, indx1, 0) m[0, 1] = np.take(thk, indx2, 0) @@ -147,14 +149,14 @@ def cascade(hk, J=7): m *= s2 # construct the grid of points - x = np.arange(0, N*(1 << J), dtype=np.float) / (1 << J) - phi = 0*x + x = np.arange(0, N * (1 << J), dtype=np.float) / (1 << J) + phi = 0 * x - psi = 0*x + psi = 0 * x # find phi0, and phi1 lam, v = eig(m[0, 0]) - ind = np.argmin(np.absolute(lam-1)) + ind = np.argmin(np.absolute(lam - 1)) # a dictionary with a binary representation of the # evaluation points x < 1 -- i.e. position is 0.xxxx v = np.real(v[:, ind]) @@ -169,29 +171,29 @@ def cascade(hk, J=7): bitdic['1'] = np.dot(m[0, 1], bitdic['0']) step = 1 << J phi[::step] = bitdic['0'] - phi[(1<<(J-1))::step] = bitdic['1'] + phi[(1 << (J - 1))::step] = bitdic['1'] psi[::step] = np.dot(m[1, 0], bitdic['0']) - psi[(1 << (J-1))::step] = np.dot(m[1, 1], bitdic['0']) + psi[(1 << (J - 1))::step] = np.dot(m[1, 1], bitdic['0']) # descend down the levels inserting more and more values # into bitdic -- store the values in the correct location once we # have computed them -- stored in the dictionary # for quicker use later. prevkeys = ['1'] - for level in range(2, J+1): + for level in range(2, J + 1): newkeys = ['%d%s' % (xx, yy) for xx in [0, 1] for yy in prevkeys] - fac = 1 << (J-level) + fac = 1 << (J - level) for key in newkeys: # convert key to number num = 0 for pos in range(level): if key[pos] == '1': - num += (1 << (level-1-pos)) + num += (1 << (level - 1 - pos)) pastphi = bitdic[key[1:]] ii = int(key[0]) temp = np.dot(m[0, ii], pastphi) bitdic[key] = temp - phi[num*fac::step] = temp - psi[num*fac::step] = np.dot(m[1, ii], pastphi) + phi[num * fac::step] = temp + psi[num * fac::step] = np.dot(m[1, ii], pastphi) prevkeys = newkeys return x, phi, psi @@ -235,12 +237,12 @@ def morlet(M, w=5.0, s=1.0, complete=True): by f = 2*s*w*r / M where r is the sampling rate. """ - x = linspace(-s*2*pi, s*2*pi, M) - output = exp(1j*w*x) + x = linspace(-s * 2 * pi, s * 2 * pi, M) + output = exp(1j * w * x) if complete: - output -= exp(-0.5*(w**2)) + output -= exp(-0.5 * (w ** 2)) - output *= exp(-0.5*(x**2)) * pi**(-0.25) + output *= exp(-0.5 * (x ** 2)) * pi ** (-0.25) return output diff --git a/scipy/signal/windows.py b/scipy/signal/windows.py index d2cc3dcc5f7c..4fb3e4c4e77d 100644 --- a/scipy/signal/windows.py +++ b/scipy/signal/windows.py @@ -28,12 +28,12 @@ def triang(M, sym=True): odd = M % 2 if not sym and not odd: M = M + 1 - n = np.arange(1, int((M+1)/2)+1) + n = np.arange(1, int((M + 1) / 2) + 1) if M % 2 == 0: - w = (2*n-1.0)/M + w = (2 * n - 1.0) / M w = np.r_[w, w[::-1]] else: - w = 2*n/(M+1.0) + w = 2 * n / (M + 1.0) w = np.r_[w, w[-2::-1]] if not sym and not odd: @@ -52,11 +52,12 @@ def parzen(M, sym=True): odd = M % 2 if not sym and not odd: M = M + 1 - n = np.arange(-(M-1)/2.0, (M-1)/2.0+0.5, 1.0) - na = np.extract(n < -(M-1)/4.0, n) - nb = np.extract(abs(n) <= (M-1)/4.0, n) - wa = 2*(1-np.abs(na)/(M/2.0))**3.0 - wb = 1-6*(np.abs(nb)/(M/2.0))**2.0 + 6*(np.abs(nb)/(M/2.0))**3.0 + n = np.arange(-(M - 1) / 2.0, (M - 1) / 2.0 + 0.5, 1.0) + na = np.extract(n < -(M - 1) / 4.0, n) + nb = np.extract(abs(n) <= (M - 1) / 4.0, n) + wa = 2 * (1 - np.abs(na) / (M / 2.0)) ** 3.0 + wb = (1 - 6 * (np.abs(nb) / (M / 2.0)) ** 2.0 + + 6 * (np.abs(nb) / (M / 2.0)) ** 3.0) w = np.r_[wa, wb, wa[::-1]] if not sym and not odd: w = w[:-1] @@ -75,7 +76,7 @@ def bohman(M, sym=True): if not sym and not odd: M = M + 1 fac = np.abs(np.linspace(-1, 1, M)[1:-1]) - w = (1 - fac) * np.cos(np.pi*fac) + 1.0/np.pi*np.sin(np.pi*fac) + w = (1 - fac) * np.cos(np.pi * fac) + 1.0 / np.pi * np.sin(np.pi * fac) w = np.r_[0, w, 0] if not sym and not odd: w = w[:-1] @@ -94,7 +95,8 @@ def blackman(M, sym=True): if not sym and not odd: M = M + 1 n = np.arange(0, M) - w = 0.42-0.5*np.cos(2.0*np.pi*n/(M-1)) + 0.08*np.cos(4.0*np.pi*n/(M-1)) + w = (0.42 - 0.5 * np.cos(2.0 * np.pi * n / (M - 1)) + + 0.08 * np.cos(4.0 * np.pi * n / (M - 1))) if not sym and not odd: w = w[:-1] return w @@ -113,8 +115,9 @@ def nuttall(M, sym=True): M = M + 1 a = [0.3635819, 0.4891775, 0.1365995, 0.0106411] n = np.arange(0, M) - fac = n*2*np.pi/(M-1.0) - w = a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac) - a[3]*np.cos(3*fac) + fac = n * 2 * np.pi / (M - 1.0) + w = (a[0] - a[1] * np.cos(fac) + + a[2] * np.cos(2 * fac) - a[3] * np.cos(3 * fac)) if not sym and not odd: w = w[:-1] return w @@ -133,8 +136,9 @@ def blackmanharris(M, sym=True): M = M + 1 a = [0.35875, 0.48829, 0.14128, 0.01168] n = np.arange(0, M) - fac = n*2*np.pi/(M-1.0) - w = a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac) - a[3]*np.cos(3*fac) + fac = n * 2 * np.pi / (M - 1.0) + w = (a[0] - a[1] * np.cos(fac) + + a[2] * np.cos(2 * fac) - a[3] * np.cos(3 * fac)) if not sym and not odd: w = w[:-1] return w @@ -153,9 +157,10 @@ def flattop(M, sym=True): M = M + 1 a = [0.2156, 0.4160, 0.2781, 0.0836, 0.0069] n = np.arange(0, M) - fac = n*2*np.pi/(M-1.0) - w = a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac) - a[3]*np.cos(3*fac) + \ - a[4]*np.cos(4*fac) + fac = n * 2 * np.pi / (M - 1.0) + w = (a[0] - a[1] * np.cos(fac) + + a[2] * np.cos(2 * fac) - a[3] * np.cos(3 * fac) + + a[4] * np.cos(4 * fac)) if not sym and not odd: w = w[:-1] return w @@ -173,7 +178,8 @@ def bartlett(M, sym=True): if not sym and not odd: M = M + 1 n = np.arange(0, M) - w = np.where(np.less_equal(n, (M-1)/2.0), 2.0*n/(M-1), 2.0-2.0*n/(M-1)) + w = np.where(np.less_equal(n, (M - 1) / 2.0), + 2.0 * n / (M - 1), 2.0 - 2.0 * n / (M - 1)) if not sym and not odd: w = w[:-1] return w @@ -191,7 +197,7 @@ def hanning(M, sym=True): if not sym and not odd: M = M + 1 n = np.arange(0, M) - w = 0.5-0.5*np.cos(2.0*np.pi*n/(M-1)) + w = 0.5 - 0.5 * np.cos(2.0 * np.pi * n / (M - 1)) if not sym and not odd: w = w[:-1] return w @@ -211,8 +217,8 @@ def barthann(M, sym=True): if not sym and not odd: M = M + 1 n = np.arange(0, M) - fac = np.abs(n/(M-1.0)-0.5) - w = 0.62 - 0.48*fac + 0.38*np.cos(2*np.pi*fac) + fac = np.abs(n / (M - 1.0) - 0.5) + w = 0.62 - 0.48 * fac + 0.38 * np.cos(2 * np.pi * fac) if not sym and not odd: w = w[:-1] return w @@ -230,7 +236,7 @@ def hamming(M, sym=True): if not sym and not odd: M = M + 1 n = np.arange(0, M) - w = 0.54-0.46*np.cos(2.0*np.pi*n/(M-1)) + w = 0.54 - 0.46 * np.cos(2.0 * np.pi * n / (M - 1)) if not sym and not odd: w = w[:-1] return w @@ -248,8 +254,9 @@ def kaiser(M, beta, sym=True): if not sym and not odd: M = M + 1 n = np.arange(0, M) - alpha = (M-1)/2.0 - w = special.i0(beta * np.sqrt(1-((n-alpha)/alpha)**2.0))/special.i0(beta) + alpha = (M - 1) / 2.0 + w = (special.i0(beta * np.sqrt(1 - ((n - alpha) / alpha) ** 2.0)) / + special.i0(beta)) if not sym and not odd: w = w[:-1] return w @@ -266,9 +273,9 @@ def gaussian(M, std, sym=True): odd = M % 2 if not sym and not odd: M = M + 1 - n = np.arange(0, M) - (M-1.0)/2.0 - sig2 = 2*std*std - w = np.exp(-n**2 / sig2) + n = np.arange(0, M) - (M - 1.0) / 2.0 + sig2 = 2 * std * std + w = np.exp(-n ** 2 / sig2) if not sym and not odd: w = w[:-1] return w @@ -288,8 +295,8 @@ def general_gaussian(M, p, sig, sym=True): odd = M % 2 if not sym and not odd: M = M + 1 - n = np.arange(0, M) - (M-1.0)/2.0 - w = np.exp(-0.5*(n/sig)**(2*p)) + n = np.arange(0, M) - (M - 1.0) / 2.0 + w = np.exp(-0.5 * (n / sig) ** (2 * p)) if not sym and not odd: w = w[:-1] return w @@ -321,15 +328,15 @@ def chebwin(M, at, sym=True): # compute the parameter beta order = M - 1.0 - beta = np.cosh(1.0/order * np.arccosh(10**(np.abs(at)/20.))) - k = np.r_[0:M]*1.0 - x = beta * np.cos(np.pi*k/M) - #find the window's DFT coefficients + beta = np.cosh(1.0 / order * np.arccosh(10 ** (np.abs(at) / 20.))) + k = np.r_[0:M] * 1.0 + x = beta * np.cos(np.pi * k / M) + # Find the window's DFT coefficients # Use analytic definition of Chebyshev polynomial instead of expansion # from scipy.special. Using the expansion in scipy.special leads to errors. p = np.zeros(x.shape) p[x > 1] = np.cosh(order * np.arccosh(x[x > 1])) - p[x < -1] = (1 - 2*(order%2)) * np.cosh(order * np.arccosh(-x[x < -1])) + p[x < -1] = (1 - 2 * (order % 2)) * np.cosh(order * np.arccosh(-x[x < -1])) p[np.abs(x) <= 1] = np.cos(order * np.arccos(x[np.abs(x) <= 1])) # Appropriate IDFT and filling up @@ -340,7 +347,7 @@ def chebwin(M, at, sym=True): w = w[:n] / w[0] w = np.concatenate((w[n - 1:0:-1], w)) else: - p = p * np.exp(1.j*np.pi / M * np.r_[0:M]) + p = p * np.exp(1.j * np.pi / M * np.r_[0:M]) w = np.real(fft(p)) n = M / 2 + 1 w = w / w[1] @@ -354,7 +361,7 @@ def slepian(M, width, sym=True): """Return the M-point slepian window. """ - if (M*width > 27.38): + if (M * width > 27.38): raise ValueError("Cannot reliably obtain slepian sequences for" " M*width > 27.38.") if M < 1: @@ -365,12 +372,12 @@ def slepian(M, width, sym=True): if not sym and not odd: M = M + 1 - twoF = width/2.0 - alpha = (M-1)/2.0 + twoF = width / 2.0 + alpha = (M - 1) / 2.0 m = np.arange(0, M) - alpha n = m[:, np.newaxis] k = m[np.newaxis, :] - AF = twoF*special.sinc(twoF*(n-k)) + AF = twoF * special.sinc(twoF * (n - k)) [lam, vec] = linalg.eig(AF) ind = np.argmax(abs(lam), axis=-1) w = np.abs(vec[:, ind])