git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/python/gnuradio/gruimpl/freqz.py

   1 #!/usr/bin/env python
   2 #
   3 # Copyright 2005,2007 Free Software Foundation, Inc.
   4 #
   5 # This file is part of GNU Radio
   6 #
   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)
  10 # any later version.
  11 #
  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.
  16 #
  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.
  21 #
  22
  23 # This code lifted from various parts of www.scipy.org -eb 2005-01-24
  24
  25 # Copyright (c) 2001, 2002 Enthought, Inc.
  26 #
  27 # All rights reserved.
  28 #
  29 # Redistribution and use in source and binary forms, with or without
  30 # modification, are permitted provided that the following conditions are met:
  31 #
  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.
  40 #
  41 #
  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
  52 # DAMAGE.
  53 #
  54
  55 __all__ = ['freqz']
  56
  57 import numpy
  58 from numpy import *
  59 Num=numpy
  60
  61 def atleast_1d(*arys):
  62     """ Force a sequence of arrays to each be at least 1D.
  63
  64          Description:
  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.
  68          Arguments:
  69             *arys -- arrays to be converted to 1 or more dimensional array.
  70          Returns:
  71             input array converted to at least 1D array.
  72     """
  73     res = []
  74     for ary in arys:
  75         ary = asarray(ary)
  76         if len(ary.shape) == 0:
  77             result = numpy.array([ary[0]])
  78         else:
  79             result = ary
  80         res.append(result)
  81     if len(res) == 1:
  82         return res[0]
  83     else:
  84         return res
  85
  86
  87 def polyval(p,x):
  88     """Evaluate the polynomial p at x.  If x is a polynomial then composition.
  89
  90     Description:
  91
  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]
  94
  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
  97       returned.
  98     """
  99     p = asarray(p)
 100     if isinstance(x,poly1d):
 101         y = 0
 102     else:
 103         x = asarray(x)
 104         y = numpy.zeros(x.shape,x.typecode())
 105     for i in range(len(p)):
 106         y = x * y + p[i]
 107     return y
 108
 109 class poly1d:
 110     """A one-dimensional polynomial class.
 111
 112     p = poly1d([1,2,3]) constructs the polynomial x**2 + 2 x + 3
 113
 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.
 120
 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.
 125     """
 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]
 130             return
 131         if r:
 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')
 137         if len(c_or_r) == 0:
 138             c_or_r = numpy.array([0])
 139         self.__dict__['coeffs'] = c_or_r
 140         self.__dict__['order'] = len(c_or_r) - 1
 141
 142     def __array__(self,t=None):
 143         if t:
 144             return asarray(self.coeffs,t)
 145         else:
 146             return asarray(self.coeffs)
 147
 148     def __coerce__(self,other):
 149         return None
 150
 151     def __repr__(self):
 152         vals = repr(self.coeffs)
 153         vals = vals[6:-1]
 154         return "poly1d(%s)" % vals
 155
 156     def __len__(self):
 157         return self.order
 158
 159     def __str__(self):
 160         N = self.order
 161         thestr = "0"
 162         for k in range(len(self.coeffs)):
 163             coefstr ='%.4g' % abs(self.coeffs[k])
 164             if coefstr[-4:] == '0000':
 165                 coefstr = coefstr[:-5]
 166             power = (N-k)
 167             if power == 0:
 168                 if coefstr != '0':
 169                     newstr = '%s' % (coefstr,)
 170                 else:
 171                     if k == 0:
 172                         newstr = '0'
 173                     else:
 174                         newstr = ''
 175             elif power == 1:
 176                 if coefstr == '0':
 177                     newstr = ''
 178                 elif coefstr == '1':
 179                     newstr = 'x'
 180                 else:
 181                     newstr = '%s x' % (coefstr,)
 182             else:
 183                 if coefstr == '0':
 184                     newstr = ''
 185                 elif coefstr == '1':
 186                     newstr = 'x**%d' % (power,)
 187                 else:
 188                     newstr = '%s x**%d' % (coefstr, power)
 189
 190             if k > 0:
 191                 if newstr != '':
 192                     if self.coeffs[k] < 0:
 193                         thestr = "%s - %s" % (thestr, newstr)
 194                     else:
 195                         thestr = "%s + %s" % (thestr, newstr)
 196             elif (k == 0) and (newstr != '') and (self.coeffs[k] < 0):
 197                 thestr = "-%s" % (newstr,)
 198             else:
 199                 thestr = newstr
 200         return _raise_power(thestr)
 201
 202
 203     def __call__(self, val):
 204         return polyval(self.coeffs, val)
 205
 206     def __mul__(self, other):
 207         if isscalar(other):
 208             return poly1d(self.coeffs * other)
 209         else:
 210             other = poly1d(other)
 211             return poly1d(polymul(self.coeffs, other.coeffs))
 212
 213     def __rmul__(self, other):
 214         if isscalar(other):
 215             return poly1d(other * self.coeffs)
 216         else:
 217             other = poly1d(other)
 218             return poly1d(polymul(self.coeffs, other.coeffs))
 219
 220     def __add__(self, other):
 221         other = poly1d(other)
 222         return poly1d(polyadd(self.coeffs, other.coeffs))
 223
 224     def __radd__(self, other):
 225         other = poly1d(other)
 226         return poly1d(polyadd(self.coeffs, other.coeffs))
 227
 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."
 231         res = [1]
 232         for k in range(val):
 233             res = polymul(self.coeffs, res)
 234         return poly1d(res)
 235
 236     def __sub__(self, other):
 237         other = poly1d(other)
 238         return poly1d(polysub(self.coeffs, other.coeffs))
 239
 240     def __rsub__(self, other):
 241         other = poly1d(other)
 242         return poly1d(polysub(other.coeffs, self.coeffs))
 243
 244     def __div__(self, other):
 245         if isscalar(other):
 246             return poly1d(self.coeffs/other)
 247         else:
 248             other = poly1d(other)
 249             return map(poly1d,polydiv(self.coeffs, other.coeffs))
 250
 251     def __rdiv__(self, other):
 252         if isscalar(other):
 253             return poly1d(other/self.coeffs)
 254         else:
 255             other = poly1d(other)
 256             return map(poly1d,polydiv(other.coeffs, self.coeffs))
 257
 258     def __setattr__(self, key, val):
 259         raise ValueError, "Attributes cannot be changed this way."
 260
 261     def __getattr__(self, key):
 262         if key in ['r','roots']:
 263             return roots(self.coeffs)
 264         elif key in ['c','coef','coefficients']:
 265             return self.coeffs
 266         elif key in ['o']:
 267             return self.order
 268         else:
 269             return self.__dict__[key]
 270
 271     def __getitem__(self, val):
 272         ind = self.order - val
 273         if val > self.order:
 274             return 0
 275         if val < 0:
 276             return 0
 277         return self.coeffs[ind]
 278
 279     def __setitem__(self, key, val):
 280         ind = self.order - key
 281         if key < 0:
 282             raise ValueError, "Does not support negative powers."
 283         if key > self.order:
 284             zr = numpy.zeros(key-self.order,self.coeffs.typecode())
 285             self.__dict__['coeffs'] = numpy.concatenate((zr,self.coeffs))
 286             self.__dict__['order'] = key
 287             ind = 0
 288         self.__dict__['coeffs'][ind] = val
 289         return
 290
 291     def integ(self, m=1, k=0):
 292         return poly1d(polyint(self.coeffs,m=m,k=k))
 293
 294     def deriv(self, m=1):
 295         return poly1d(polyder(self.coeffs,m=m))
 296
 297 def freqz(b, a, worN=None, whole=0, plot=None):
 298     """Compute frequency response of a digital filter.
 299
 300     Description:
 301
 302        Given the numerator (b) and denominator (a) of a digital filter compute
 303        its frequency response.
 304
 305                   jw               -jw            -jmw
 306            jw  B(e)    b[0] + b[1]e + .... + b[m]e
 307         H(e) = ---- = ------------------------------------
 308                   jw               -jw            -jnw
 309                A(e)    a[0] + a[2]e + .... + a[n]e
 310
 311     Inputs:
 312
 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
 319                 to 2*pi.
 320
 321     Outputs: (h,w)
 322
 323        h -- The frequency response.
 324        w -- The frequencies at which h was computed.
 325     """
 326     b, a = map(atleast_1d, (b,a))
 327     if whole:
 328         lastpoint = 2*pi
 329     else:
 330         lastpoint = pi
 331     if worN is None:
 332         N = 512
 333         w = Num.arange(0,lastpoint,lastpoint/N)
 334     elif isinstance(worN, types.IntType):
 335         N = worN
 336         w = Num.arange(0,lastpoint,lastpoint/N)
 337     else:
 338         w = worN
 339     w = atleast_1d(w)
 340     zm1 = exp(-1j*w)
 341     h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
 342     # if not plot is None:
 343     #    plot(w, h)
 344     return h, w