problem with numpy

Jeff Boyd boyd at cpsc.ucalgary.ca
Fri Jan 5 12:29:49 EST 2001


I have a program that uses Numerical Python that is behaving oddly. A
program that demonstrates is below.  The assignment to y[i] in the
function "filter" (an IIR filter) does not appear to happen so the
result is always an array of zeros.  Questions:

1. is this a bug or am I doing something wrong?
2.  if it is a bug, will an upgrade fix it? (I am using ver 15.2 right
now)



# bug.py - a bug, I think, in NumPy

import Numeric
from math import *

# apply an IIR filter (H(z) = N(z)/D(z)) to an input signal, x
def filter( N, D, x ) :
    y = Numeric.zeros( x.shape )        # output array
    n = x.shape[0]                      # number of samples
    N = N / D[0]                        # normalize filter coefficients
    D = D / D[0]
    on = N.shape[0] - 1                 # orders of N and D polynomials
    od = D.shape[0] - 1
    N = N[-1::-1]                       # reverse order of coefficients
    D = D[-1::-1]
    for i in range(1,n) :
        #print i, Numeric.sum( x[i-on:i+1] * N ),
        #print Numeric.sum( y[i-od:i  ] * D[:-1] )

        #############################
        # this assignment fails
        y[i] = Numeric.sum( x[i-on:i+1] * N ) - Numeric.sum( y[i-od:i  ]
* D[:-1] )
        #############################

    return y

# compute filter coefficients
tau = 5
w1 = 2 * tan(pi/tau)
Nlp = Numeric.array( [1.0, 1.0] )
Dlp = Numeric.array( [1.0 + 2/w1, 1.0 - 2/w1] )

# apply the filter to a step function
x = Numeric.zeros( 30 )
x[15:] = 1
y = filter( Nlp, Dlp, x )

# print the result
for i in range(30) :
    print i, y[i]



--
------------------------------------------------------------------------------
Jeffrey E.  Boyd
Department of Computer Science               E-mail: boyd at cpsc.ucalgary.ca
University of Calgary                        Phone: (403) 220 6038
2500 University Drive NW                     Fax: (403) 284 4707
Calgary, Alberta, Canada T2N 1N4             URL: http://www.cpsc.ucalgary.ca/
------------------------------------------------------------------------------






More information about the Python-list mailing list