frequency analysis without numpy

debug domelectric at gmail.com
Tue Jan 20 23:09:33 CET 2009


Hi-

I've been using python now for about 2 months for plugin development
within Maya (a commercial 3d application). I'm currently in the
process of writing a sound analysis plugin for maya and have completed
a good portion of it including the ability to retrieve the amplitude
at certain intervals from a wav file.

The next part of the project however consists of analysing the wav
file and outputting the amplitude at certain frequencies. Due to the
need of distributing this to many computers after it has been finished
I'm reluctant to use an external python module for FFT'ing. Accuracy
and speed aren't really issues as just a generalisation of the
amplitude of low, medium and high frequencies is required.

Do you either know of any generic FFT functions written in python that
could just be inserted into my processing class (open source
licensed).

So far i've managed to put together a chunk of code but I'm not sure
its returning the right values, any ideas?

import wave
import sys
import array
import struct
from cmath import pi, exp


def nextpow2(i):
   n = 2
   while n < i: n = n * 2
   return n

def bitrev(x):
   N, x = len(x), x[:]
   if N != nextpow2(N): raise ValueError, 'N is not power of 2'
   for i in range(N):
      k, b, a = 0, N>>1, 1
      while b >= a:
         if b & i: k = k | a
         if a & i: k = k | b
         b, a = b>>1, a<<1
   if i < k:      # important not to swap back
      x[i], x[k] = x[k], x[i]
   return x



def fft(x, sign=-1):
   N, W = len(x), []
   for i in range(N):      # exp(-j...) is default
      W.append(exp(sign * 2j * pi * i / N))
   x = bitrev(x)
   m = 2
   while m <= N:
      for s in range(0, N, m):
         for i in range(m/2):
            n = i * N / m
            a, b = s + i, s + i + m/2
            x[a], x[b] = x[a] + W[n % N] * x[b], x[a] - W[n % N] * x
[b]
      m = m * 2
   return x



def ifft(X):
    N, x = len(X), fft(X, sign=1)   # e^{j2\pi/N}
    for i in range(N):
        x[i] = x[i] / float(N)
    return x




fp = wave.open("hankuncom.wav","rb")
sample_rate = fp.getframerate()
total_num_samps = fp.getnframes()
fft_length = 4096
num_fft = (total_num_samps / fft_length ) - 2


temp = []

for i in range(num_fft):
   tempb = fp.readframes(fft_length)
   moo = struct.unpack("%uh"%(fft_length),tempb)
   temp.append(moo)


listed = list(temp[0])

freq_pwr  = fft(listed)

print(freq_pwr)

fp.close()

Quick warning this code snippet does run slowly
Thanks!

The wav file referenced in the script is here: www.reality-debug.co.uk/hankuncom.zip



More information about the Python-list mailing list