summaryrefslogtreecommitdiff
path: root/gnuradio-core/src/python/gnuradio/gruimpl/freqz.py
diff options
context:
space:
mode:
Diffstat (limited to 'gnuradio-core/src/python/gnuradio/gruimpl/freqz.py')
-rw-r--r--gnuradio-core/src/python/gnuradio/gruimpl/freqz.py344
1 files changed, 0 insertions, 344 deletions
diff --git a/gnuradio-core/src/python/gnuradio/gruimpl/freqz.py b/gnuradio-core/src/python/gnuradio/gruimpl/freqz.py
deleted file mode 100644
index 60dca64a58..0000000000
--- a/gnuradio-core/src/python/gnuradio/gruimpl/freqz.py
+++ /dev/null
@@ -1,344 +0,0 @@
-#!/usr/bin/env python
-#
-# Copyright 2005,2007 Free Software Foundation, Inc.
-#
-# This file is part of GNU Radio
-#
-# GNU Radio is free software; you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation; either version 3, or (at your option)
-# any later version.
-#
-# GNU Radio is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with GNU Radio; see the file COPYING. If not, write to
-# the Free Software Foundation, Inc., 51 Franklin Street,
-# Boston, MA 02110-1301, USA.
-#
-
-# This code lifted from various parts of www.scipy.org -eb 2005-01-24
-
-# Copyright (c) 2001, 2002 Enthought, Inc.
-#
-# All rights reserved.
-#
-# Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions are met:
-#
-# a. Redistributions of source code must retain the above copyright notice,
-# this list of conditions and the following disclaimer.
-# b. Redistributions in binary form must reproduce the above copyright
-# notice, this list of conditions and the following disclaimer in the
-# documentation and/or other materials provided with the distribution.
-# c. Neither the name of the Enthought nor the names of its contributors
-# may be used to endorse or promote products derived from this software
-# without specific prior written permission.
-#
-#
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-# ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
-# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
-# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
-# DAMAGE.
-#
-
-__all__ = ['freqz']
-
-import numpy
-from numpy import *
-Num=numpy
-
-def atleast_1d(*arys):
- """ Force a sequence of arrays to each be at least 1D.
-
- Description:
- Force an array to be at least 1D. If an array is 0D, the
- array is converted to a single row of values. Otherwise,
- the array is unaltered.
- Arguments:
- *arys -- arrays to be converted to 1 or more dimensional array.
- Returns:
- input array converted to at least 1D array.
- """
- res = []
- for ary in arys:
- ary = asarray(ary)
- if len(ary.shape) == 0:
- result = numpy.array([ary[0]])
- else:
- result = ary
- res.append(result)
- if len(res) == 1:
- return res[0]
- else:
- return res
-
-
-def polyval(p,x):
- """Evaluate the polynomial p at x. If x is a polynomial then composition.
-
- Description:
-
- If p is of length N, this function returns the value:
- p[0]*(x**N-1) + p[1]*(x**N-2) + ... + p[N-2]*x + p[N-1]
-
- x can be a sequence and p(x) will be returned for all elements of x.
- or x can be another polynomial and the composite polynomial p(x) will be
- returned.
- """
- p = asarray(p)
- if isinstance(x,poly1d):
- y = 0
- else:
- x = asarray(x)
- y = numpy.zeros(x.shape,x.typecode())
- for i in range(len(p)):
- y = x * y + p[i]
- return y
-
-class poly1d:
- """A one-dimensional polynomial class.
-
- p = poly1d([1,2,3]) constructs the polynomial x**2 + 2 x + 3
-
- p(0.5) evaluates the polynomial at the location
- p.r is a list of roots
- p.c is the coefficient array [1,2,3]
- p.order is the polynomial order (after leading zeros in p.c are removed)
- p[k] is the coefficient on the kth power of x (backwards from
- sequencing the coefficient array.
-
- polynomials can be added, substracted, multplied and divided (returns
- quotient and remainder).
- asarray(p) will also give the coefficient array, so polynomials can
- be used in all functions that accept arrays.
- """
- def __init__(self, c_or_r, r=0):
- if isinstance(c_or_r,poly1d):
- for key in c_or_r.__dict__.keys():
- self.__dict__[key] = c_or_r.__dict__[key]
- return
- if r:
- c_or_r = poly(c_or_r)
- c_or_r = atleast_1d(c_or_r)
- if len(c_or_r.shape) > 1:
- raise ValueError, "Polynomial must be 1d only."
- c_or_r = trim_zeros(c_or_r, trim='f')
- if len(c_or_r) == 0:
- c_or_r = numpy.array([0])
- self.__dict__['coeffs'] = c_or_r
- self.__dict__['order'] = len(c_or_r) - 1
-
- def __array__(self,t=None):
- if t:
- return asarray(self.coeffs,t)
- else:
- return asarray(self.coeffs)
-
- def __coerce__(self,other):
- return None
-
- def __repr__(self):
- vals = repr(self.coeffs)
- vals = vals[6:-1]
- return "poly1d(%s)" % vals
-
- def __len__(self):
- return self.order
-
- def __str__(self):
- N = self.order
- thestr = "0"
- for k in range(len(self.coeffs)):
- coefstr ='%.4g' % abs(self.coeffs[k])
- if coefstr[-4:] == '0000':
- coefstr = coefstr[:-5]
- power = (N-k)
- if power == 0:
- if coefstr != '0':
- newstr = '%s' % (coefstr,)
- else:
- if k == 0:
- newstr = '0'
- else:
- newstr = ''
- elif power == 1:
- if coefstr == '0':
- newstr = ''
- elif coefstr == '1':
- newstr = 'x'
- else:
- newstr = '%s x' % (coefstr,)
- else:
- if coefstr == '0':
- newstr = ''
- elif coefstr == '1':
- newstr = 'x**%d' % (power,)
- else:
- newstr = '%s x**%d' % (coefstr, power)
-
- if k > 0:
- if newstr != '':
- if self.coeffs[k] < 0:
- thestr = "%s - %s" % (thestr, newstr)
- else:
- thestr = "%s + %s" % (thestr, newstr)
- elif (k == 0) and (newstr != '') and (self.coeffs[k] < 0):
- thestr = "-%s" % (newstr,)
- else:
- thestr = newstr
- return _raise_power(thestr)
-
-
- def __call__(self, val):
- return polyval(self.coeffs, val)
-
- def __mul__(self, other):
- if isscalar(other):
- return poly1d(self.coeffs * other)
- else:
- other = poly1d(other)
- return poly1d(polymul(self.coeffs, other.coeffs))
-
- def __rmul__(self, other):
- if isscalar(other):
- return poly1d(other * self.coeffs)
- else:
- other = poly1d(other)
- return poly1d(polymul(self.coeffs, other.coeffs))
-
- def __add__(self, other):
- other = poly1d(other)
- return poly1d(polyadd(self.coeffs, other.coeffs))
-
- def __radd__(self, other):
- other = poly1d(other)
- return poly1d(polyadd(self.coeffs, other.coeffs))
-
- def __pow__(self, val):
- if not isscalar(val) or int(val) != val or val < 0:
- raise ValueError, "Power to non-negative integers only."
- res = [1]
- for k in range(val):
- res = polymul(self.coeffs, res)
- return poly1d(res)
-
- def __sub__(self, other):
- other = poly1d(other)
- return poly1d(polysub(self.coeffs, other.coeffs))
-
- def __rsub__(self, other):
- other = poly1d(other)
- return poly1d(polysub(other.coeffs, self.coeffs))
-
- def __div__(self, other):
- if isscalar(other):
- return poly1d(self.coeffs/other)
- else:
- other = poly1d(other)
- return map(poly1d,polydiv(self.coeffs, other.coeffs))
-
- def __rdiv__(self, other):
- if isscalar(other):
- return poly1d(other/self.coeffs)
- else:
- other = poly1d(other)
- return map(poly1d,polydiv(other.coeffs, self.coeffs))
-
- def __setattr__(self, key, val):
- raise ValueError, "Attributes cannot be changed this way."
-
- def __getattr__(self, key):
- if key in ['r','roots']:
- return roots(self.coeffs)
- elif key in ['c','coef','coefficients']:
- return self.coeffs
- elif key in ['o']:
- return self.order
- else:
- return self.__dict__[key]
-
- def __getitem__(self, val):
- ind = self.order - val
- if val > self.order:
- return 0
- if val < 0:
- return 0
- return self.coeffs[ind]
-
- def __setitem__(self, key, val):
- ind = self.order - key
- if key < 0:
- raise ValueError, "Does not support negative powers."
- if key > self.order:
- zr = numpy.zeros(key-self.order,self.coeffs.typecode())
- self.__dict__['coeffs'] = numpy.concatenate((zr,self.coeffs))
- self.__dict__['order'] = key
- ind = 0
- self.__dict__['coeffs'][ind] = val
- return
-
- def integ(self, m=1, k=0):
- return poly1d(polyint(self.coeffs,m=m,k=k))
-
- def deriv(self, m=1):
- return poly1d(polyder(self.coeffs,m=m))
-
-def freqz(b, a, worN=None, whole=0, plot=None):
- """Compute frequency response of a digital filter.
-
- Description:
-
- Given the numerator (b) and denominator (a) of a digital filter compute
- its frequency response.
-
- jw -jw -jmw
- jw B(e) b[0] + b[1]e + .... + b[m]e
- H(e) = ---- = ------------------------------------
- jw -jw -jnw
- A(e) a[0] + a[2]e + .... + a[n]e
-
- Inputs:
-
- b, a --- the numerator and denominator of a linear filter.
- worN --- If None, then compute at 512 frequencies around the unit circle.
- If a single integer, the compute at that many frequencies.
- Otherwise, compute the response at frequencies given in worN
- whole -- Normally, frequencies are computed from 0 to pi (upper-half of
- unit-circle. If whole is non-zero compute frequencies from 0
- to 2*pi.
-
- Outputs: (h,w)
-
- h -- The frequency response.
- w -- The frequencies at which h was computed.
- """
- b, a = map(atleast_1d, (b,a))
- if whole:
- lastpoint = 2*pi
- else:
- lastpoint = pi
- if worN is None:
- N = 512
- w = Num.arange(0,lastpoint,lastpoint/N)
- elif isinstance(worN, types.IntType):
- N = worN
- w = Num.arange(0,lastpoint,lastpoint/N)
- else:
- w = worN
- w = atleast_1d(w)
- zm1 = exp(-1j*w)
- h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
- # if not plot is None:
- # plot(w, h)
- return h, w