3 # Copyright 2005,2007 Free Software Foundation, Inc.
5 # This file is part of GNU Radio
7 # GNU Radio is free software; you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation; either version 3, or (at your option)
12 # GNU Radio is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
17 # You should have received a copy of the GNU General Public License
18 # along with GNU Radio; see the file COPYING. If not, write to
19 # the Free Software Foundation, Inc., 51 Franklin Street,
20 # Boston, MA 02110-1301, USA.
23 # This code lifted from various parts of www.scipy.org -eb 2005-01-24
25 # Copyright (c) 2001, 2002 Enthought, Inc.
27 # All rights reserved.
29 # Redistribution and use in source and binary forms, with or without
30 # modification, are permitted provided that the following conditions are met:
32 # a. Redistributions of source code must retain the above copyright notice,
33 # this list of conditions and the following disclaimer.
34 # b. Redistributions in binary form must reproduce the above copyright
35 # notice, this list of conditions and the following disclaimer in the
36 # documentation and/or other materials provided with the distribution.
37 # c. Neither the name of the Enthought nor the names of its contributors
38 # may be used to endorse or promote products derived from this software
39 # without specific prior written permission.
42 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
43 # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
44 # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
45 # ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
46 # ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
47 # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
48 # SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
49 # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
50 # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
51 # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
61 def atleast_1d(*arys):
62 """ Force a sequence of arrays to each be at least 1D.
65 Force an array to be at least 1D. If an array is 0D, the
66 array is converted to a single row of values. Otherwise,
67 the array is unaltered.
69 *arys -- arrays to be converted to 1 or more dimensional array.
71 input array converted to at least 1D array.
76 if len(ary.shape) == 0:
77 result = numpy.array([ary[0]])
88 """Evaluate the polynomial p at x. If x is a polynomial then composition.
92 If p is of length N, this function returns the value:
93 p[0]*(x**N-1) + p[1]*(x**N-2) + ... + p[N-2]*x + p[N-1]
95 x can be a sequence and p(x) will be returned for all elements of x.
96 or x can be another polynomial and the composite polynomial p(x) will be
100 if isinstance(x,poly1d):
104 y = numpy.zeros(x.shape,x.typecode())
105 for i in range(len(p)):
110 """A one-dimensional polynomial class.
112 p = poly1d([1,2,3]) constructs the polynomial x**2 + 2 x + 3
114 p(0.5) evaluates the polynomial at the location
115 p.r is a list of roots
116 p.c is the coefficient array [1,2,3]
117 p.order is the polynomial order (after leading zeros in p.c are removed)
118 p[k] is the coefficient on the kth power of x (backwards from
119 sequencing the coefficient array.
121 polynomials can be added, substracted, multplied and divided (returns
122 quotient and remainder).
123 asarray(p) will also give the coefficient array, so polynomials can
124 be used in all functions that accept arrays.
126 def __init__(self, c_or_r, r=0):
127 if isinstance(c_or_r,poly1d):
128 for key in c_or_r.__dict__.keys():
129 self.__dict__[key] = c_or_r.__dict__[key]
132 c_or_r = poly(c_or_r)
133 c_or_r = atleast_1d(c_or_r)
134 if len(c_or_r.shape) > 1:
135 raise ValueError, "Polynomial must be 1d only."
136 c_or_r = trim_zeros(c_or_r, trim='f')
138 c_or_r = numpy.array([0])
139 self.__dict__['coeffs'] = c_or_r
140 self.__dict__['order'] = len(c_or_r) - 1
142 def __array__(self,t=None):
144 return asarray(self.coeffs,t)
146 return asarray(self.coeffs)
148 def __coerce__(self,other):
152 vals = repr(self.coeffs)
154 return "poly1d(%s)" % vals
162 for k in range(len(self.coeffs)):
163 coefstr ='%.4g' % abs(self.coeffs[k])
164 if coefstr[-4:] == '0000':
165 coefstr = coefstr[:-5]
169 newstr = '%s' % (coefstr,)
181 newstr = '%s x' % (coefstr,)
186 newstr = 'x**%d' % (power,)
188 newstr = '%s x**%d' % (coefstr, power)
192 if self.coeffs[k] < 0:
193 thestr = "%s - %s" % (thestr, newstr)
195 thestr = "%s + %s" % (thestr, newstr)
196 elif (k == 0) and (newstr != '') and (self.coeffs[k] < 0):
197 thestr = "-%s" % (newstr,)
200 return _raise_power(thestr)
203 def __call__(self, val):
204 return polyval(self.coeffs, val)
206 def __mul__(self, other):
208 return poly1d(self.coeffs * other)
210 other = poly1d(other)
211 return poly1d(polymul(self.coeffs, other.coeffs))
213 def __rmul__(self, other):
215 return poly1d(other * self.coeffs)
217 other = poly1d(other)
218 return poly1d(polymul(self.coeffs, other.coeffs))
220 def __add__(self, other):
221 other = poly1d(other)
222 return poly1d(polyadd(self.coeffs, other.coeffs))
224 def __radd__(self, other):
225 other = poly1d(other)
226 return poly1d(polyadd(self.coeffs, other.coeffs))
228 def __pow__(self, val):
229 if not isscalar(val) or int(val) != val or val < 0:
230 raise ValueError, "Power to non-negative integers only."
233 res = polymul(self.coeffs, res)
236 def __sub__(self, other):
237 other = poly1d(other)
238 return poly1d(polysub(self.coeffs, other.coeffs))
240 def __rsub__(self, other):
241 other = poly1d(other)
242 return poly1d(polysub(other.coeffs, self.coeffs))
244 def __div__(self, other):
246 return poly1d(self.coeffs/other)
248 other = poly1d(other)
249 return map(poly1d,polydiv(self.coeffs, other.coeffs))
251 def __rdiv__(self, other):
253 return poly1d(other/self.coeffs)
255 other = poly1d(other)
256 return map(poly1d,polydiv(other.coeffs, self.coeffs))
258 def __setattr__(self, key, val):
259 raise ValueError, "Attributes cannot be changed this way."
261 def __getattr__(self, key):
262 if key in ['r','roots']:
263 return roots(self.coeffs)
264 elif key in ['c','coef','coefficients']:
269 return self.__dict__[key]
271 def __getitem__(self, val):
272 ind = self.order - val
277 return self.coeffs[ind]
279 def __setitem__(self, key, val):
280 ind = self.order - key
282 raise ValueError, "Does not support negative powers."
284 zr = numpy.zeros(key-self.order,self.coeffs.typecode())
285 self.__dict__['coeffs'] = numpy.concatenate((zr,self.coeffs))
286 self.__dict__['order'] = key
288 self.__dict__['coeffs'][ind] = val
291 def integ(self, m=1, k=0):
292 return poly1d(polyint(self.coeffs,m=m,k=k))
294 def deriv(self, m=1):
295 return poly1d(polyder(self.coeffs,m=m))
297 def freqz(b, a, worN=None, whole=0, plot=None):
298 """Compute frequency response of a digital filter.
302 Given the numerator (b) and denominator (a) of a digital filter compute
303 its frequency response.
306 jw B(e) b[0] + b[1]e + .... + b[m]e
307 H(e) = ---- = ------------------------------------
309 A(e) a[0] + a[2]e + .... + a[n]e
313 b, a --- the numerator and denominator of a linear filter.
314 worN --- If None, then compute at 512 frequencies around the unit circle.
315 If a single integer, the compute at that many frequencies.
316 Otherwise, compute the response at frequencies given in worN
317 whole -- Normally, frequencies are computed from 0 to pi (upper-half of
318 unit-circle. If whole is non-zero compute frequencies from 0
323 h -- The frequency response.
324 w -- The frequencies at which h was computed.
326 b, a = map(atleast_1d, (b,a))
333 w = Num.arange(0,lastpoint,lastpoint/N)
334 elif isinstance(worN, types.IntType):
336 w = Num.arange(0,lastpoint,lastpoint/N)
341 h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
342 # if not plot is None: