[SciPy-user] Help with array functions
eric jones
eric at enthought.com
Tue May 4 12:31:52 EDT 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:
> Yikes. Your right!
>
> 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
More information about the SciPy-User
mailing list