[python]DFT(discrete fourier transform) and FFT
"""DFT and FFT"""
import math
def iexp(n):
return complex(math.cos(n), math.sin(n))
def is_pow2(n):
return False if n == 0 else (n == 1 or is_pow2(n >> 1))
def dft(xs):
"naive dft"
n = len(xs)
return [sum((xs[k] * iexp(-2 * math.pi * i * k / n) for k in range(n)))
for i in range(n)]
def dftinv(xs):
"naive dft"
n = len(xs)
return [sum((xs[k] * iexp(2 * math.pi * i * k / n) for k in range(n))) / n
for i in range(n)]
def fft_(xs, n, start=0, stride=1):
"cooley-turkey fft"
if n == 1: return [xs[start]]
hn, sd = n // 2, stride * 2
rs = fft_(xs, hn, start, sd) + fft_(xs, hn, start + stride, sd)
for i in range(hn):
e = iexp(-2 * math.pi * i / n)
rs[i], rs[i + hn] = rs[i] + e * rs[i + hn], rs[i] - e * rs[i + hn]
pass
return rs
def fft(xs):
assert is_pow2(len(xs))
return fft_(xs, len(xs))
def fftinv_(xs, n, start=0, stride=1):
"cooley-turkey fft"
if n == 1: return [xs[start]]
hn, sd = n // 2, stride * 2
rs = fftinv_(xs, hn, start, sd) + fftinv_(xs, hn, start + stride, sd)
for i in range(hn):
e = iexp(2 * math.pi * i / n)
rs[i], rs[i + hn] = rs[i] + e * rs[i + hn], rs[i] - e * rs[i + hn]
pass
return rs
def fftinv(xs):
assert is_pow2(len(xs))
n = len(xs)
return [v / n for v in fftinv_(xs, n)]
if __name__ == "__main__":
wave = [0, 1, 2, 3, 3, 2, 1, 0]
dfreq = dft(wave)
ffreq = fft(wave)
dwave = dftinv(dfreq)
fwave= fftinv(ffreq)
print(dfreq)
print(ffreq)
print([v.real for v in dwave])
print([v.real for v in fwave])
pass
import sys
from numpy import NaN, Inf, arange, isscalar, asarray, array
def peakdet(v, delta, x = None):
"""
Converted from MATLAB script at http://billauer.co.il/peakdet.html
Returns two arrays
function [maxtab, mintab]=peakdet(v, delta, x)
%PEAKDET Detect peaks in a vector
% [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
% maxima and minima ("peaks") in the vector V.
% MAXTAB and MINTAB consists of two columns. Column 1
% contains indices in V, and column 2 the found values.
%
% With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
% in MAXTAB and MINTAB are replaced with the corresponding
% X-values.
%
% A point is considered a maximum peak if it has the maximal
% value, and was preceded (to the left) by a value lower by
% DELTA.
% Eli Billauer, 3.4.05 (Explicitly not copyrighted).
% This function is released to the public domain; Any use is allowed.
"""
maxtab = []
mintab = []
if x is None:
x = arange(len(v))
v = asarray(v)
if len(v) != len(x):
sys.exit('Input vectors v and x must have same length')
if not isscalar(delta):
sys.exit('Input argument delta must be a scalar')
if delta <= 0:
sys.exit('Input argument delta must be positive')
mn, mx = Inf, -Inf
mnpos, mxpos = NaN, NaN
lookformax = True
for i in arange(len(v)):
this = v[i]
if this > mx:
mx = this
mxpos = x[i]
if this < mn:
mn = this
mnpos = x[i]
if lookformax:
if this < mx-delta:
maxtab.append((mxpos, mx))
mn = this
mnpos = x[i]
lookformax = False
else:
if this > mn+delta:
mintab.append((mnpos, mn))
mx = this
mxpos = x[i]
lookformax = True
return array(maxtab), array(mintab)
if __name__=="__main__":
from matplotlib.pyplot import plot, scatter, show
series = [0,0,0,2,0,0,0,-2,0,0,0,2,0,0,0,-2,0]
maxtab, mintab = peakdet(series,.3)
plot(series)
scatter(array(maxtab)[:,0], array(maxtab)[:,1], color='blue')
scatter(array(mintab)[:,0], array(mintab)[:,1], color='red')
show()