NumPy Question: Audio Frequency Spectrum / Signal Analysis - pe_audioLevels.py [0/1]

[Fernando Pérez]

Perry Kivolowitz wrote:

> [This followup was posted to comp.lang.python and a copy was sent to the
> [cited author.]
> 
> I'm having a little touble writing a graphical view of an audio
> spectrum. I'm using NumPy, and Useful Things - a Python programmable
> plug-in for Adobe After Effects.
> 
> I am totally inexperienced with respect to signal processing, so please
> excuse my (possibly) bogus use of terminology.
> 
> It seems to work, with the exception that sounds close to silence
> (highly negative DB) are putting out enormous amplitudes. So much so, I
> have to put in a base DB check - squelching anything more negative to 0.
> 
> This is the code which converts from samples to frequency data:
> 
> m = 2.0 / float(numSamples)
> halfSamples = numSamples / 2
> samplePeriod = 1.0 / float(sampleRate) * numSamples
> firstHarmonic = 1.0 / samplePeriod
> fftData = fft(samples)

> realTerms = fftData[1 : halfSamples + 1]
  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Unnecessary. Use FFT.real_fft() which takes into account the periodicity of 
the fft of a real signal and only returns the first half. Just remember to 
invert with FFT.inverse_real_fft()


> then I do:
> 
> frequencyAmplitudes = [ 0 ] * halfSamples
> for i in xrange(halfSamples):
> frequencyAmplitudes[i] = convertToDecibels(m * abs(realTerms[i]))

This is probably taking forever and  a day. Write it as a single Numeric 
array operation. As much as possible, try to avoid explicit python loops over 
numeric arrays like the plague. We're talking factors of 100-1000 in 
performance difference here.

Something like (untested) should do:

realTerms = FFT.real_fft(samples)
frequencyAmplitudes = toDB(m*Numeric.absolute(realTerms))

> to gather up data from a specific band I finish with (to show the min DB
> test):
> 
> outputData[i] = 1.0 - (abs(bandAccumulator / bandCounter) - minDB)
> / (maxDB - minDB)
> if outputData[i] > 1: outputData[i] = 1
> if outputData[i] < 0: outputData[i] = 0

Again, no looping. Use Numeric.clip() for this.

> where convertToDecibels is:
> 
> amp2db = 20.0 / log(10.0)
> 
> def convertToDecibels(a):
> if a < 1.0e-30: return -600
> else: return amp2db * log(a)

This can be written as a Numeric operation with no looping. Just mask out the 
values below your clipping threshold with where() or clip() so your log 
doesn't blow up.

> When the min DB setting is played with, I get great results. Without it
> being set right, I get blown out in the lower frequencies.
> 
> If this isn't enough, I have enclosed the entire source. If you're an
> After Effects user, I can supply a free license to Useful Things so you
> can actually see what's going on (and code your own After Effects plug-
> ins in Python!) in exchange for your help.
> 
> Thanks

cheers,

f