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

[eric jones]

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