# [SciPy-user] Help with array functions

eric jones eric at enthought.com
Tue May 4 11:31:52 CDT 2004

So, try two...

I don't know of a pure Python solution either.  My next option is
usually playing around in weave.  Here is a weave snippet that
implements the algorithm in C.  It runs, but just like the last one is
untested.  Given the track record here, a decent test is advisable... :-)

Here is the output from my cursory check to see if the code would work:

C:\temp>wv.py
weave results:  [ 0. 0.8 0.96 0.992 0.9984  0.99968  0.999936
0.9999872  0.99999744  1.]
python results: [ 0. 0.8 0.96 0.992 0.9984  0.99968  0.999936
0.9999872  0.99999744  1.]
C time: 0.00707567096418
Python time: 0.316835262617
weave vs. Python speedup: 44.7781226997

This is using gcc as the compiler in weave and 2048 samples in a
signal.  If you drop to 256, the improvement is cut in half.

Now the code:

# test.py
from scipy import zeros, ones, Float64, arange
from scipy import stats
import weave

def python_algorithm(gain, 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, 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', '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
weave_algorithm(gain, complexout, attack, peak, level, samples)
print 'weave results:', gain[:10]

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

# speed tests
import time

N = 100

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

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

eric jones wrote:

>
> eric
>
> Perry Greenfield wrote:
>
>> Eric Jones wrote:
>>
>>
>>> Bob.Cowdery at CGI-Europe.com wrote:
>>>
>>> The general ideas of uisng Numeric in Python is to avoid loops, even
>>> at the expense of sometimes doing extra computations.  In this case,
>>> the signal processing calculations can be done for every entry in
>>> the array (except the first).  The conditions checked for in the
>>> "if" statements can then be filtered out in subsequent array
>>> processing steps.
>>>
>>>
>>> Here is some (totally untested) code that should give you some ideas:
>>>
>>>    gain = self.m_procGain
>>>    abs_complexout = abs(self.m_complexout)
>>>
>>>    # Do singal processing calculation for all elements in array
>>>    gain[1:] = ( (gain[:-1] * (1 - attack)) +
>>>                 ((attack * peak) / abs_complexout[:-1]) )
>>>
>>>
>>
>> Does this do what is needed? I think the point he was making
>> was that gain[i] depends on gain[i-1] which in turn depends
>> on gain[i-2] and so on. The expression above doesn't capture that
>> since from his code he doesn't show any initial value for gain
>> (implicitly the gain[0] must exist). So using an array expression
>> which uses a shifted gain array in its input won't work.
>> This is the sort of problem that doesn't fit the array model
>> very well unless there are some primitive array operations
>> that happen to do the iterative thing needed (sometimes
>> the accumulate or convolve functions can be used for this
>> thing, but I can't see any obvious existing function. Seems
>> like a candidate for a C extension but maybe I'm missing something.
>>
>> Perry
>>
>>
>>
>> _______________________________________________
>> SciPy-user mailing list
>> SciPy-user at scipy.net
>> http://www.scipy.net/mailman/listinfo/scipy-user
>>
>>
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user