***[Possible UCE]*** [SciPy-user] Help with array functions
eric jones
eric at enthought.com
Tue May 4 15:25:34 CDT 2004
Hey Travis,
The Fortran version didn't work for me. It looks like g77 is detected
on my machine, and I get compile errors in the code. I also have an f90
compiler on my machine that wasn't picked up. Is there a way to switch
which compiler is used through command line?
eric
Travis Oliphant wrote:
> 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
>>
>>
>
> _______________________________________________
> 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