# Python frequency analysis without numpy

Started by January 20, 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):
moo = struct.unpack("%uh"%(fft_length),tempb)
temp.append(moo)

listed = list(temp[0])

freq_pwr  = fft(listed)

print(freq_pwr)

fp.close()

Thanks!


debug <dom@reality-debug.co.uk> wrote:

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

Have you thought about using boost.python?

Steve

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

Even though you don't want to use numpy in the distribution use it to
verify you functions.  You can create a unit test that will compare
the results of your function to those of fftw.  You can build in
tolerated error into your unit test.

Example of unit test here, http://www.myhdl.org/doku.php/users:cfelton:projects:recursivefft

>debug <dom@reality-debug.co.uk> wrote:
>
>>Do you either know of any generic FFT functions written in python that
>>could just be inserted into my processing class (open source licensed).
>
>Have you thought about using boost.python?
>
>Steve
>

The boost library contains FFT routines?

If you don't want to use numpy you can allways right your own
wrapper for kissfft, djbfft, or FFTW (or any other fft library you can
find).

If you want a pure python implementation (I wouldn't know why):