[SciPy-dev] problem with signal.butter

John Hunter jdhunter at ace.bsd.uchicago.edu
Thu Aug 4 07:56:20 CDT 2005


>>>>> "Fernando" == Fernando Perez <Fernando.Perez at colorado.edu> writes:

    Fernando> By all means do.  0.3.2 is quite out of date, so there's
    Fernando> a good chance CVS will take care of things.

That did the trick.  Because I need this fix for software that I
distribute to windoze users, I still have some interest in seeing if I
can backport the fix the 0.3.2, since this is the current binary
release and thus I suspect it will also be the release in the new
enthought python.  I'm posting the diff betwween scipy 0.3.3_309.4626
and 0.3.2 in the signal module, in case the person(s) who made those
changes will have some insight into what is contributing to the fix.
If not, it's no big deal because I can always try and build scipy CVS
on windows (arggg) or just tell folks they can't use this feature.  It
would be much easier just to say, "drop this file in
site-packages/scipy/signal" if it's just a one or two line fix.

Thanks!
JDH

Only in scipy/Lib/signal/: CVS
Only in scipy/Lib/signal/: docs
diff -cr SciPy_complete-0.3.2/Lib/signal/filter_design.py scipy/Lib/signal/filter_design.py
*** SciPy_complete-0.3.2/Lib/signal/filter_design.py	2004-02-26 02:55:12.000000000 -0600
--- scipy/Lib/signal/filter_design.py	2004-11-22 14:25:00.000000000 -0600
***************
*** 5,10 ****
--- 5,11 ----
  from scipy_base.fastumath import *
  from scipy_base import atleast_1d, poly, polyval, roots, imag, real, asarray,\
       allclose, Float, resize, pi, concatenate, absolute, logspace, c_
+ from scipy_base import mintypecode
  from scipy import comb, special, optimize, linalg
  import string, types
  
***************
*** 228,236 ****
      a,b = map(atleast_1d,(a,b))
      D = len(a) - 1
      N = len(b) - 1
!     artype = b.typecode()
!     if artype not in ['F','D','f','d']:
!         artype = Num.Float
      ma = max([N,D])
      Np = N + ma
      Dp = D + ma
--- 229,235 ----
      a,b = map(atleast_1d,(a,b))
      D = len(a) - 1
      N = len(b) - 1
!     artype = mintypecode((a,b))
      ma = max([N,D])
      Np = N + ma
      Dp = D + ma
***************
*** 261,269 ****
      a,b = map(atleast_1d,(a,b))
      D = len(a) - 1
      N = len(b) - 1
!     artype = b.typecode()
!     if artype not in ['F','D','f','d']:
!         artype = Num.Float
      M = max([N,D])
      Np = M + M
      Dp = M + M
--- 260,266 ----
      a,b = map(atleast_1d,(a,b))
      D = len(a) - 1
      N = len(b) - 1
!     artype = mintypecode((a,b))
      M = max([N,D])
      Np = M + M
      Dp = M + M
diff -cr SciPy_complete-0.3.2/Lib/signal/firfilter.c scipy/Lib/signal/firfilter.c
*** SciPy_complete-0.3.2/Lib/signal/firfilter.c	2003-07-08 20:52:20.000000000 -0500
--- scipy/Lib/signal/firfilter.c	2005-05-31 01:51:32.000000000 -0500
***************
*** 112,118 ****
    if ((sum = calloc(type_size,2))==NULL) return -3; /* No memory */
    value = sum + type_size;
  
!   if (outsize == FULL) {Os[0] = Ns[0]+Nwin[0]-1; Os[1] = Ns[0]+Nwin[1]-1;}
    else if (outsize == SAME) {Os[0] = Ns[0]; Os[1] = Ns[1];}
    else if (outsize == VALID) {Os[0] = Ns[0]-Nwin[0]+1; Os[1] = Ns[1]-Nwin[1]+1;}
    else return -1;  /* Invalid output flag */  
--- 112,118 ----
    if ((sum = calloc(type_size,2))==NULL) return -3; /* No memory */
    value = sum + type_size;
  
!   if (outsize == FULL) {Os[0] = Ns[0]+Nwin[0]-1; Os[1] = Ns[1]+Nwin[1]-1;}
    else if (outsize == SAME) {Os[0] = Ns[0]; Os[1] = Ns[1];}
    else if (outsize == VALID) {Os[0] = Ns[0]-Nwin[0]+1; Os[1] = Ns[1]-Nwin[1]+1;}
    else return -1;  /* Invalid output flag */  
***************
*** 120,125 ****
--- 120,128 ----
    if ((boundary != PAD) && (boundary != REFLECT) && (boundary != CIRCULAR)) 
      return -2;   /* Invalid boundary flag */
  
+   /* Speed this up by not doing any if statements in the for loop.  Need 3*3*2=18 different
+      loops executed for different conditions */
+ 
    for (m=0; m < Os[0]; m++) {
      /* Reposition index into input image based on requested output size */
      if (outsize == FULL) new_m = convolve ? m : (m-Nwin[0]+1);
diff -cr SciPy_complete-0.3.2/Lib/signal/info_signal.py scipy/Lib/signal/info_signal.py
*** SciPy_complete-0.3.2/Lib/signal/info_signal.py	2004-02-26 01:15:10.000000000 -0600
--- scipy/Lib/signal/info_signal.py	2005-05-31 01:51:32.000000000 -0500
***************
*** 6,15 ****
--- 6,17 ----
   
      convolve      --  N-dimensional convolution.
      correlate     --  N-dimensional correlation.
+     fftconvolve   --  N-dimensional convolution using the FFT.    
      convolve2d    --  2-dimensional convolution (more options).
      correlate2d   --  2-dimensional correlation (more options).
      sepfir2d      --  Convolve with a 2-D separable FIR filter.
  
+ 
   B-splines:
      
      bspline       --  B-spline basis function of order n.
***************
*** 45,52 ****
      freqs         --- Analog filter frequency response
      freqz         --- Digital filter frequency response
  
! 
!  matlab-style IIR filter design:
   
      butter (buttord)  -- Butterworth
      cheby1 (cheb1ord) -- Chebyshev Type I
--- 47,53 ----
      freqs         --- Analog filter frequency response
      freqz         --- Digital filter frequency response
  
!  Matlab-style IIR filter design:
   
      butter (buttord)  -- Butterworth
      cheby1 (cheb1ord) -- Chebyshev Type I
***************
*** 68,74 ****
      tf2ss -- transfer function to state-space.
      ss2tf -- state-pace to transfer function.
      zpk2ss -- zero-pole-gain to state-space.
!     ss2zpk -- state-space to pole-zero-gain.    
  """
  
  postpone_import = 1
--- 69,88 ----
      tf2ss -- transfer function to state-space.
      ss2tf -- state-pace to transfer function.
      zpk2ss -- zero-pole-gain to state-space.
!     ss2zpk -- state-space to pole-zero-gain.
! 
!  Waveforms:
! 
!     sawtooth -- Periodic sawtooth
!     square -- Square wave
!     gausspulse -- Gaussian modulated sinusoid
!     chirp -- Frequency swept cosine signal
! 
!  Wavelets:
! 
!     daub -- return low-pass filter for daubechies wavelets
!     qmf  -- return quadrature mirror filter from low-pass
!     cascade -- compute scaling function and wavelet from coefficients
  """
  
  postpone_import = 1
diff -cr SciPy_complete-0.3.2/Lib/signal/__init__.py scipy/Lib/signal/__init__.py
*** SciPy_complete-0.3.2/Lib/signal/__init__.py	2004-04-13 16:23:51.000000000 -0500
--- scipy/Lib/signal/__init__.py	2005-05-31 01:51:32.000000000 -0500
***************
*** 10,13 ****
--- 10,14 ----
  from filter_design import *
  from ltisys import *
  from signaltools import *
+ from wavelets import *
  
diff -cr SciPy_complete-0.3.2/Lib/signal/ltisys.py scipy/Lib/signal/ltisys.py
*** SciPy_complete-0.3.2/Lib/signal/ltisys.py	2003-12-13 07:44:51.000000000 -0600
--- scipy/Lib/signal/ltisys.py	2005-02-24 04:42:24.000000000 -0600
***************
*** 223,230 ****
              self.__dict__['zeros'], self.__dict__['poles'], \
                                      self.__dict__['gain'] = ss2zpk(*args)
              self.__dict__['num'], self.__dict__['den'] = ss2tf(*args)
!             self.inputs = B.shape[-1]
!             self.outputs = C.shape[0]
          else:
              raise ValueError, "Needs 2, 3, or 4 arguments."
  
--- 223,230 ----
              self.__dict__['zeros'], self.__dict__['poles'], \
                                      self.__dict__['gain'] = ss2zpk(*args)
              self.__dict__['num'], self.__dict__['den'] = ss2tf(*args)
!             self.inputs = self.B.shape[-1]
!             self.outputs = self.C.shape[0]
          else:
              raise ValueError, "Needs 2, 3, or 4 arguments."
  
diff -cr SciPy_complete-0.3.2/Lib/signal/setup_signal.py scipy/Lib/signal/setup_signal.py
*** SciPy_complete-0.3.2/Lib/signal/setup_signal.py	2004-02-06 02:17:52.000000000 -0600
--- scipy/Lib/signal/setup_signal.py	2005-03-24 04:22:16.000000000 -0600
***************
*** 5,10 ****
--- 5,11 ----
  from scipy_distutils.misc_util import get_path, default_config_dict, dot_join
  
  def configuration(parent_package='',parent_path=None):
+     from scipy_distutils.system_info import get_info
      package = 'signal'
      local_path = get_path(__name__,parent_path)    
      config = default_config_dict(package, parent_package)
diff -cr SciPy_complete-0.3.2/Lib/signal/signaltools.py scipy/Lib/signal/signaltools.py
*** SciPy_complete-0.3.2/Lib/signal/signaltools.py	2004-06-02 17:57:55.000000000 -0500
--- scipy/Lib/signal/signaltools.py	2005-01-17 18:53:50.000000000 -0600
***************
*** 10,17 ****
  import types
  import scipy
  from scipy.stats import mean
! import Numeric
! from Numeric import array, arange, where, sqrt, rank, zeros, NewAxis, argmax
  from scipy_base.fastumath import *
  
  _modedict = {'valid':0, 'same':1, 'full':2}
--- 10,17 ----
  import types
  import scipy
  from scipy.stats import mean
! import scipy_base as Numeric
! from scipy_base import array, arange, where, sqrt, rank, zeros, NewAxis, argmax, product
  from scipy_base.fastumath import *
  
  _modedict = {'valid':0, 'same':1, 'full':2}
***************
*** 66,72 ****
      kernel = asarray(in2)
      if rank(volume) == rank(kernel) == 0:
          return volume*kernel
!     if (Numeric.product(kernel.shape) > Numeric.product(volume.shape)):
          temp = kernel
          kernel = volume
          volume = temp
--- 66,72 ----
      kernel = asarray(in2)
      if rank(volume) == rank(kernel) == 0:
          return volume*kernel
!     if (product(kernel.shape) > product(volume.shape)):
          temp = kernel
          kernel = volume
          volume = temp
***************
*** 76,81 ****
--- 76,117 ----
  
      return sigtools._correlateND(volume, kernel, val)
  
+ def _centered(arr, newsize):
+     # Return the center newsize portion of the array.
+     newsize = asarray(newsize)
+     currsize = array(arr.shape)
+     startind = (currsize - newsize) / 2
+     endind = startind + newsize
+     myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
+     return arr[tuple(myslice)]
+         
+ def fftconvolve(in1, in2, mode="full"):
+     """Convolve two N-dimensional arrays using FFT. SEE convolve
+     """
+     s1 = array(in1.shape)
+     s2 = array(in2.shape)
+     if (s1.typecode() in ['D','F']) or (s2.typecode() in ['D', 'F']):
+         cmplx=1
+     else: cmplx=0
+     size = s1+s2-1
+     IN1 = fftn(in1,size)
+     IN1 *= fftn(in2,size)
+     ret = ifftn(IN1)
+     del IN1
+     if not cmplx:
+         ret = real(ret)
+     if mode == "full":
+         return ret
+     elif mode == "same":
+         if product(s1) > product(s2):
+             osize = s1
+         else:
+             osize = s2
+         return _centered(ret,osize)
+     elif mode == "valid":
+         return _centered(ret,abs(s2-s1)+1)
+             
+ 
  def convolve(in1, in2, mode='full'):
      """Convolve two N-dimensional arrays.
  
Only in scipy/Lib/signal/tests: CVS
Only in scipy/Lib/signal/: wavelets.py




More information about the Scipy-dev mailing list