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

Travis Oliphant oliphant at ee.byu.edu
Tue May 4 13:47:43 CDT 2004

```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
>
>

```