***[Possible UCE]*** [SciPy-user] Help with array functions

[Travis Oliphant]

Bob.Cowdery at CGI-Europe.com wrote:

> Hi all,
>  
> Can anyone help me convert this loop to an array function. I can't 
> have loops in my signal processing as they take far too much 
> processing power. This is a fragment from a speech processor using 
> Ferkins formula. Essentially the gain array (m_procGain) which is 2048 
> long (m_sampl_sz) needs its nth element processing into it's (n-1)th 
> element. This is eventually used to multiply the samples in 
> m_complexout, the first 2048 elements of this are used. Attack is a 
> scaler currently 0.4.
>
>                 peak = 
> abs(self.m_complexout[argmax(abs(self.m_complexout[:self.m_sampl_sz]))])
>                 i = arange(1,self.m_sampl_sz)
>                 for x in i:
>                     if abs(self.m_complexout[x-1]) > 0:
>                        self.m_procGain[x] = ( (self.m_procGain[x-1] * 
> (1 - attack)) + ((attack * peak) / abs(self.m_complexout[x-1])) )        
>                        if self.m_procGain[x] > level:
>                            self.m_procGain[x] = level
>                     else:
>                        self.m_procGain[x] = 1
>

Bob,

You have here something similar to an IIR filter.  The command 
signal.lfilter is usually used in situations like this.

If I read this code correctly then define y[n] = self.m_procGain[n]  
with a = (1-attack)

w[n] = (attack*peak) / abs(self.m_complexout[n])

and you want to basically compute

y[n] = a y[n-1] + w[n-1]  

as long as w[n-1] is not infinite because of zero-valued 
abs(self.m_complexout[n-1])

The problem with using signal.lfilter is you have a non-linear filter 
due to the if statements.  You could try to modify the algorithm a bit 
to maintain a linear filter, or you will have to use weave (or f2py). 

Just to expose you to the many ways to accomplish the problem.  Here is 
another test using f2py



******** Results *********************
weave results: [ 0.      0.8     0.96    0.992   0.9984  1.      1.      
1.      1.   
       1.    ]
python results: [ 0.      0.8     0.96    0.992   0.9984  1.      
1.      1.      1.   
       1.    ]
f2py results: [ 0.      0.8     0.96    0.992   0.9984  1.      1.      
1.      1.   
       1.    ]
C time: 0.08
Python time: 6.13
f2py time: 0.05
weave vs. Python speedup: 76.625
f2py vs. Python speedup: 122.6


********* Code below ***********



# test.py

fort_code = r"""

       subroutine gain_compute(gain, abs_complexout, attack, peak,
     $        level, samples)

       integer samples
       double precision attack, peak, level
       double precision gain(samples), abs_complexout(samples)

       double precision val
       integer i

       do i = 2, samples
           if (abs_complexout(i-1).gt.0.0) then
               val = gain(i-1)*(1.0-attack) +
     $                 (attack*peak) / abs_complexout(i-1)
               if (val.gt.level) then
                   val = level
               end if
           else
              val = 1.0
           end if
           gain(i) = val
       end do

       return
       end subroutine gain_compute
"""
     

from scipy import zeros, ones, Float64, arange
from scipy import stats
import scipy_distutils
import weave
import f2py2e

try:
    import forttest
except ImportError:
    f2py2e.compile(fort_code, modulename='forttest')
    import forttest

def python_algorithm(gain, abs_complexout, attack, peak, level, samples):
   # gain is changed *inplace*
   i = arange(1,samples)
   for x in i:
       if abs_complexout[x-1] > 0:
          gain[x] = ( (gain[x-1] * (1.0 - attack)) + ((attack * peak) / 
abs_complexout[x-1]) )
          if gain[x] > level:
              gain[x] = level
       else:
          gain[x] = 1

def weave_algorithm(gain, abs_complexout, attack, peak, level, samples):
   # gain is changed *inplace*
   code = """
          for (int x=1; x < samples; x++)
          {
              if (abs_complexout[x-1] > 0.0)
              {
                 gain[x] = ( (gain[x-1] * (1.0 - attack)) +
                             ((attack * peak) / abs_complexout[x-1]) 
);                        if (gain[x] > level)
                     gain[x] = level;
              }
              else
                 gain[x] = 1;
          }
          """
   weave.inline(code,
                ['gain', 'abs_complexout', 'samples', 'attack', 'peak', 
'level'],
                compiler='gcc')
 
# some default parameters...
samples = 2048
gain = zeros(samples,typecode=Float64)
complexout = stats.choice([-1.0,0.0,1.0],size=samples)

attack = .8
peak = 1.0
level = 1.0


# equivalence tests
ac = abs(complexout)
weave_algorithm(gain, ac, attack, peak, level, samples)
print 'weave results:', gain[:10]

gain = zeros(samples,typecode=Float64)
python_algorithm(gain, ac, attack, peak, level, samples)
print 'python results:', gain[:10]

gain = zeros(samples,typecode=Float64)
forttest.gain_compute(gain, ac, attack, peak, level, samples)
print 'f2py results:', gain[:10]


# speed tests
import time

N = 1000

t1 = time.clock()
for it in range(N):
   weave_algorithm(gain, complexout, attack, peak, level, samples)
t2 = time.clock()
ct = t2 - t1
print 'C time:', ct

t1 = time.clock()
for it in range(N):
   python_algorithm(gain, complexout, attack, peak, level, samples)
t2 = time.clock()
pt = t2-t1
print 'Python time:', pt

f2py_algorithm = forttest.gain_compute
t1 = time.clock()
for it in range(N):
   f2py_algorithm(gain, complexout, attack, peak, level, samples)
t2 = time.clock()
ft = t2-t1
print 'f2py time:', ft


ratio = pt/ct
ratio2 = pt/ft
print 'weave vs. Python speedup:', ratio
print 'f2py vs. Python speedup:', ratio2












> Thanks for any assistance.
>  
> Bob
>  
>  
> Bob Cowdery
> CGI Senior Technical Architect
> +44(0)1438 791517
> Mobile: +44(0)7771 532138
> bob.cowdery at cgi-europe.com <mailto:bob.cowdery at cgi-europe.com>
>  
>  
>  
>
> **** Confidentiality Notice **** Proprietary/Confidential
> Information belonging to CGI Group Inc. and its affiliates
> may be contained in this message. If you are not a recipient
> indicated or intended in this message (or responsible for
> delivery of this message to such person), or you think for
> any reason that this message may have been addressed to you
> in error, you may not use or copy or deliver this message
> to anyone else.  In such case, you should destroy this
> message and are asked to notify the sender by reply email.
>
>------------------------------------------------------------------------
>
>_______________________________________________
>SciPy-user mailing list
>SciPy-user at scipy.net
>http://www.scipy.net/mailman/listinfo/scipy-user
>  
>