Re: frequency analysis without numpy
On Jan 21, 12:13 am, sturlamolden sturlamol...@yahoo.no wrote: Apart from that, an FFT in pure python is going to be atrociously slow for anything but the shortest signals. I cannot imagine why you want to do this. Just to elaborate on this: The whole purpose of using FFT is speed. That pretty much excludes the use of Python. If you don't care about speed, you could just as well compute the DFT directly. The FFT is just a O(n lon n) algorithm for computing the DFT. Here is a DFT with O(N**2) behavior: from math import sin, cos, pi def real_dft(x): ''' DFT for a real valued sequence x ''' r = [] N = len(x) M = N//2 + 1 if N%2 else N//2 for n in range(M): s = 0j for k in range(N): tmp = 2*pi*k*n/N s += x[k] * (cos(tmp) - 1j*sin(tmp)) r.append(s) return r S.M. -- http://mail.python.org/mailman/listinfo/python-list
Re: frequency analysis without numpy
On Jan 20, 3:13 pm, sturlamolden sturlamol...@yahoo.no wrote: Consider using Thompson's multitaper method, autoregression (maximum entropy), or Welch method for your frequency estimates. Blackman- Tuckey is also a possibility, but I see no reason to prefer that to Welch. Multitaper and AR tends to be the better options though, but they are not as fast as Welch' method. There is a fortran program named rainbow that has some dozen of these methods to compare with. Rainbow is located at http://www.digitalCalculus.com/demo/rainbow.html. I recommend Burg's method or AutoCorrelation but its data dependent for best choice. Phil -- http://mail.python.org/mailman/listinfo/python-list
Re: frequency analysis without numpy
On Jan 21, 11:31 am, sturlamolden sturlamol...@yahoo.no wrote: On Jan 21, 12:13 am, sturlamolden sturlamol...@yahoo.no wrote: Apart from that, an FFT in pure python is going to be atrociously slow for anything but the shortest signals. I cannot imagine why you want to do this. Just to elaborate on this: The whole purpose of using FFT is speed. That pretty much excludes the use of Python. If you don't care about speed, you could just as well compute the DFT directly. The FFT is just a O(n lon n) algorithm for computing the DFT. Here is a DFT with O(N**2) behavior: from math import sin, cos, pi def real_dft(x): ''' DFT for a real valued sequence x ''' r = [] N = len(x) M = N//2 + 1 if N%2 else N//2 for n in range(M): s = 0j for k in range(N): tmp = 2*pi*k*n/N s += x[k] * (cos(tmp) - 1j*sin(tmp)) r.append(s) return r S.M. Thanks for the quick reply, so what do I pass the real_dft function (obviously a list) but do I pass it the unpacked binary data? Dom -- http://mail.python.org/mailman/listinfo/python-list
frequency analysis without numpy
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, N1, 1
while b = a:
if b i: k = k | a
if a i: k = k | b
b, a = b1, a1
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
--
http://mail.python.org/mailman/listinfo/python-list
Re: frequency analysis without numpy
On Jan 20, 11:09 pm, debug domelect...@gmail.com wrote: So far i've managed to put together a chunk of code but I'm not sure its returning the right values, any ideas? Don't use the periodogram for frequency analysis; it is not a good estimate of the power spectrum. According to the Wiener-Khintchin theorem, the power spectrum is the Fourier transform of the autocorrelation, not the Fourier transform of the signal. You don't want to just FFT the signal. FFTs are used for spectral estimation, but not to obtain the periodogram. The periodogrgam is full of noise and spectral leakage. And even worse, the variance of the periodogram is constant, not inversely proportional to the number of samples. This is not how you want a spectrum estimator to behave. Consider using Thompson's multitaper method, autoregression (maximum entropy), or Welch method for your frequency estimates. Blackman- Tuckey is also a possibility, but I see no reason to prefer that to Welch. Multitaper and AR tends to be the better options though, but they are not as fast as Welch' method. You probably also want to use real FFTs instead of complex FFTs, unless you are interested in estimating the negative frequencies as well. For common spectral analysis, real FFTs are almost always what you want (complex FFTs can be used, but the incur redundant work). Apart from that, an FFT in pure python is going to be atrociously slow for anything but the shortest signals. I cannot imagine why you want to do this. There are many free FFT libraries around, including FFTPACK (used by NumPy) and FFTW. The easiest way to check your FFT is to compare with a verified implementation. I'd use NumPy or Matlab for that. S.M. -- http://mail.python.org/mailman/listinfo/python-list
