TEST_PLOTS = False
-def floor_pow_of_two(num) :
+def floor_pow_of_two(num):
"Round num down to the closest exact a power of two."
lnum = _numpy.log2(num)
- if int(lnum) != lnum :
+ if int(lnum) != lnum:
num = 2**_numpy.floor(lnum)
return num
-def round_pow_of_two(num) :
+def round_pow_of_two(num):
"Round num to the closest exact a power of two on a log scale."
lnum = _numpy.log2(num)
- if int(lnum) != lnum :
+ if int(lnum) != lnum:
num = 2**_numpy.round(lnum)
return num
-def ceil_pow_of_two(num) :
+def ceil_pow_of_two(num):
"Round num up to the closest exact a power of two."
lnum = _numpy.log2(num)
- if int(lnum) != lnum :
+ if int(lnum) != lnum:
num = 2**_numpy.ceil(lnum)
return num
-def _test_rfft(xs, Xs) :
+def _test_rfft(xs, Xs):
# Numpy's FFT algoritm returns
# n-1
# X[k] = SUM x[m] exp (-j 2pi km /n)
j = _numpy.complex(0,1)
n = len(xs)
Xa = []
- for k in range(n) :
+ for k in range(n):
Xa.append(sum([x*_numpy.exp(-j*2*_numpy.pi*k*m/n)
for x,m in zip(xs,range(n))]))
if k < len(Xs):
"Mismatch on Parseval's, {} != 1/{} * {}".format(
timeSum, n, freqSum))
-def _test_rfft_suite() :
+def _test_rfft_suite():
print('Test numpy rfft definition')
xs = [1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1]
_test_rfft(xs, _numpy.fft.rfft(xs))
-def unitary_rfft(data, freq=1.0) :
+def unitary_rfft(data, freq=1.0):
"""Compute the Fourier transform of real data.
Unitary (preserves power [Parseval's theorem]).
raise ValueError(
'Mismatch in spacing, {} != 1/({}*{})'.format(df, len(xs), dt))
Xa = list(Xs)
- for k in range(len(Xs)-1,1,-1) :
+ for k in range(len(Xs)-1,1,-1):
Xa.append(Xa[k])
if len(xs) != len(Xa):
raise ValueError('Length mismatch {} != {}'.format(len(xs), len(Xa)))
freqs,Xs = unitary_rfft(xs, 1.0/dt)
_test_unitary_rfft_parsevals(xs, 1.0/dt, freqs, Xs)
-def _rect(t) :
- if _numpy.abs(t) < 0.5 :
+def _rect(t):
+ if _numpy.abs(t) < 0.5:
return 1
- else :
+ else:
return 0
-def _test_unitary_rfft_rect(a=1.0, time_shift=5.0, samp_freq=25.6, samples=256) :
+def _test_unitary_rfft_rect(
+ a=1.0, time_shift=5.0, samp_freq=25.6, samples=256):
"Show fft(rect(at)) = 1/abs(a) * _numpy.sinc(f/a)"
samp_freq = _numpy.float(samp_freq)
a = _numpy.float(a)
x = _numpy.zeros((samples,), dtype=_numpy.float)
dt = 1.0/samp_freq
- for i in range(samples) :
+ for i in range(samples):
t = i*dt
x[i] = _rect(a*(t-time_shift))
freq_axis, X = unitary_rfft(x, samp_freq)
# remove the phase due to our time shift
j = _numpy.complex(0.0,1.0) # sqrt(-1)
- for i in range(len(freq_axis)) :
+ for i in range(len(freq_axis)):
f = freq_axis[i]
inverse_phase_shift = _numpy.exp(j*2.0*_numpy.pi*time_shift*f)
X[i] *= inverse_phase_shift
# so sinc(0.5) = sin(pi/2)/(pi/2) = 2/pi
if _numpy.sinc(0.5) != 2.0/_numpy.pi:
raise ValueError('abnormal sinc()')
- for i in range(len(freq_axis)) :
+ for i in range(len(freq_axis)):
f = freq_axis[i]
expected[i] = 1.0/_numpy.abs(a) * _numpy.sinc(f/a)
- if TEST_PLOTS :
+ if TEST_PLOTS:
figure = _pyplot.figure()
time_axes = figure.add_subplot(2, 1, 1)
time_axes.plot(_numpy.arange(0, dt*samples, dt), x)
freq_axes.plot(freq_axis, expected, 'b-')
freq_axes.set_title('freq series')
-def _test_unitary_rfft_rect_suite() :
+def _test_unitary_rfft_rect_suite():
print('Test unitary FFTs on variously shaped rectangular functions')
_test_unitary_rfft_rect(a=0.5)
_test_unitary_rfft_rect(a=2.0)
_test_unitary_rfft_rect(a=0.7, samp_freq=50, samples=512)
_test_unitary_rfft_rect(a=3.0, samp_freq=60, samples=1024)
-def _gaussian(a, t) :
+def _gaussian(a, t):
return _numpy.exp(-a * t**2)
-def _test_unitary_rfft_gaussian(a=1.0, time_shift=5.0, samp_freq=25.6, samples=256) :
+def _test_unitary_rfft_gaussian(
+ a=1.0, time_shift=5.0, samp_freq=25.6, samples=256):
"Show fft(rect(at)) = 1/abs(a) * sinc(f/a)"
samp_freq = _numpy.float(samp_freq)
a = _numpy.float(a)
x = _numpy.zeros((samples,), dtype=_numpy.float)
dt = 1.0/samp_freq
- for i in range(samples) :
+ for i in range(samples):
t = i*dt
x[i] = _gaussian(a, (t-time_shift))
freq_axis, X = unitary_rfft(x, samp_freq)
# remove the phase due to our time shift
j = _numpy.complex(0.0,1.0) # sqrt(-1)
- for i in range(len(freq_axis)) :
+ for i in range(len(freq_axis)):
f = freq_axis[i]
inverse_phase_shift = _numpy.exp(j*2.0*_numpy.pi*time_shift*f)
X[i] *= inverse_phase_shift
expected = _numpy.zeros((len(freq_axis),), dtype=_numpy.float)
- for i in range(len(freq_axis)) :
+ for i in range(len(freq_axis)):
f = freq_axis[i]
# see Wikipedia, or do the integral yourself.
expected[i] = _numpy.sqrt(_numpy.pi/a) * _gaussian(
1.0/a, _numpy.pi*f)
- if TEST_PLOTS :
+ if TEST_PLOTS:
figure = _pyplot.figure()
time_axes = figure.add_subplot(2, 1, 1)
time_axes.plot(_numpy.arange(0, dt*samples, dt), x)
freq_axes.plot(freq_axis, expected, 'b-')
freq_axes.set_title('freq series')
-def _test_unitary_rfft_gaussian_suite() :
+def _test_unitary_rfft_gaussian_suite():
print("Test unitary FFTs on variously shaped gaussian functions")
_test_unitary_rfft_gaussian(a=0.5)
_test_unitary_rfft_gaussian(a=2.0)
-def power_spectrum(data, freq=1.0) :
+def power_spectrum(data, freq=1.0):
"""Compute the power spectrum of DATA taken with a sampling frequency FREQ.
DATA must be real (not complex).
power = (trans * trans.conj()).real # We want the square of the amplitude.
return (freq_axis, power)
-def unitary_power_spectrum(data, freq=1.0) :
+def unitary_power_spectrum(data, freq=1.0):
freq_axis,power = power_spectrum(data, freq)
- # One sided power spectral density, so 2|H(f)|**2 (see NR 2nd edition 12.0.14, p498)
+ # One sided power spectral density, so 2|H(f)|**2
+ # (see NR 2nd edition 12.0.14, p498)
#
# numpy normalizes with 1/N on the inverse transform ifft,
# so we should normalize the freq-space representation with 1/sqrt(N).
- # But we're using the rfft, where N points are like N/2 complex points, so 1/sqrt(N/2)
+ # But we're using the rfft, where N points are like N/2 complex points,
+ # so 1/sqrt(N/2)
# So the power gets normalized by that twice and we have 2/N
#
# On top of this, the FFT assumes a sampling freq of 1 per second,
# and we want to preserve area under our curves.
# If our total time T = len(data)/freq is smaller than 1,
- # our df_real = freq/len(data) is bigger that the FFT expects (dt_fft = 1/len(data)),
+ # our df_real = freq/len(data) is bigger that the FFT expects
+ # (dt_fft = 1/len(data)),
# and we need to scale the powers down to conserve area.
# df_fft * F_fft(f) = df_real *F_real(f)
# F_real = F_fft(f) * (1/len)/(freq/len) = F_fft(f)*freq
return (freq_axis, power)
-def _test_unitary_power_spectrum_sin(sin_freq=10, samp_freq=512, samples=1024) :
+def _test_unitary_power_spectrum_sin(sin_freq=10, samp_freq=512, samples=1024):
x = _numpy.zeros((samples,), dtype=_numpy.float)
samp_freq = _numpy.float(samp_freq)
- for i in range(samples) :
+ for i in range(samples):
x[i] = _numpy.sin(2.0 * _numpy.pi * (i/samp_freq) * sin_freq)
freq_axis, power = unitary_power_spectrum(x, samp_freq)
imax = _numpy.argmax(power)
# P(f0) = 0.5/df
# where f0 is the sin's frequency
#
- # or :
+ # or:
# FFT of sin(2*pi*t*f0) gives
# X(f) = 0.5 i (delta(f-f0) - delta(f-f0)),
# (area under x(t) = 0, area under X(f) = 0)
sin_freq, expected[i], freq_axis[imax], power[imax]))
Pexp = 0
P = 0
- for i in range(len(freq_axis)) :
+ for i in range(len(freq_axis)):
Pexp += expected[i] *df
P += power[i] * df
print('The total power should be {} ({})'.format(Pexp, P))
- if TEST_PLOTS :
+ if TEST_PLOTS:
figure = _pyplot.figure()
time_axes = figure.add_subplot(2, 1, 1)
time_axes.plot(
freq_axes.set_title(
'{} samples of sin at {} Hz'.format(samples, sin_freq))
-def _test_unitary_power_spectrum_sin_suite() :
+def _test_unitary_power_spectrum_sin_suite():
print('Test unitary power spectrums on variously shaped sin functions')
_test_unitary_power_spectrum_sin(sin_freq=5, samp_freq=512, samples=1024)
_test_unitary_power_spectrum_sin(sin_freq=5, samp_freq=512, samples=2048)
_test_unitary_power_spectrum_sin(sin_freq=5, samp_freq=512, samples=4098)
_test_unitary_power_spectrum_sin(sin_freq=7, samp_freq=512, samples=1024)
_test_unitary_power_spectrum_sin(sin_freq=5, samp_freq=1024, samples=2048)
- # finally, with some irrational numbers, to check that I'm not getting lucky
+ # with some irrational numbers, to check that I'm not getting lucky
_test_unitary_power_spectrum_sin(
sin_freq=_numpy.pi, samp_freq=100*_numpy.exp(1), samples=1024)
# test with non-integer number of periods
_test_unitary_power_spectrum_sin(sin_freq=5, samp_freq=512, samples=256)
-def _test_unitary_power_spectrum_delta(amp=1, samp_freq=1, samples=256) :
+def _test_unitary_power_spectrum_delta(amp=1, samp_freq=1, samples=256):
x = _numpy.zeros((samples,), dtype=_numpy.float)
samp_freq = _numpy.float(samp_freq)
x[0] = amp
print('The power should be flat at y = {} ({})'.format(
expected_amp, power[0]))
- if TEST_PLOTS :
+ if TEST_PLOTS:
figure = _pyplot.figure()
time_axes = figure.add_subplot(2, 1, 1)
time_axes.plot(
freq_axes.plot(freq_axis, expected, 'b-')
freq_axes.set_title('{} samples of delta amp {}'.format(samples, amp))
-def _test_unitary_power_spectrum_delta_suite() :
+def _test_unitary_power_spectrum_delta_suite():
print('Test unitary power spectrums on various delta functions')
_test_unitary_power_spectrum_delta(amp=1, samp_freq=1.0, samples=1024)
_test_unitary_power_spectrum_delta(amp=1, samp_freq=1.0, samples=2048)
- _test_unitary_power_spectrum_delta(amp=1, samp_freq=0.5, samples=2048)# expected = 2*computed
- _test_unitary_power_spectrum_delta(amp=1, samp_freq=2.0, samples=2048)# expected = 0.5*computed
+ # expected = 2*computed
+ _test_unitary_power_spectrum_delta(amp=1, samp_freq=0.5, samples=2048)
+ # expected = 0.5*computed
+ _test_unitary_power_spectrum_delta(amp=1, samp_freq=2.0, samples=2048)
_test_unitary_power_spectrum_delta(amp=3, samp_freq=1.0, samples=1024)
_test_unitary_power_spectrum_delta(
amp=_numpy.pi, samp_freq=_numpy.exp(1), samples=1024)
-def _gaussian2(area, mean, std, t) :
+def _gaussian2(area, mean, std, t):
"Integral over all time = area (i.e. normalized for area=1)"
return area/(std*_numpy.sqrt(2.0*_numpy.pi)) * _numpy.exp(
-0.5*((t-mean)/std)**2)
-def _test_unitary_power_spectrum_gaussian(area=2.5, mean=5, std=1, samp_freq=10.24 ,samples=512) : #1024
+def _test_unitary_power_spectrum_gaussian(
+ area=2.5, mean=5, std=1, samp_freq=10.24 ,samples=512):
x = _numpy.zeros((samples,), dtype=_numpy.float)
mean = _numpy.float(mean)
- for i in range(samples) :
+ for i in range(samples):
t = i/_numpy.float(samp_freq)
x[i] = _gaussian2(area, mean, std, t)
freq_axis, power = unitary_power_spectrum(x, samp_freq)
# generate the predicted curve
# by comparing our _gaussian2() form to _gaussian(),
# we see that the Fourier transform of x(t) has parameters:
- # std' = 1/(2 pi std) (references declaring std' = 1/std are converting to angular frequency, not frequency like we are)
+ # std' = 1/(2 pi std) (references declaring std' = 1/std are
+ # converting to angular frequency,
+ # not frequency like we are)
# area' = area/[std sqrt(2*pi)] (plugging into FT of _gaussian() above)
- # mean' = 0 (changing the mean in the time-domain just changes the phase in the freq-domain)
+ # mean' = 0 (changing the mean in the time-domain just
+ # changes the phase in the freq-domain)
# So our power spectral density per unit time is given by
# P(f) = 2 |X(f)|**2 / T
# Where
expected = _numpy.zeros((len(freq_axis),), dtype=_numpy.float)
# 1/total_time ( = freq_axis[1]-freq_axis[0] = freq_axis[1])
df = _numpy.float(samp_freq)/samples
- for i in range(len(freq_axis)) :
+ for i in range(len(freq_axis)):
f = i*df
gaus = _gaussian2(area, mean, std, f)
expected[i] = 2.0 * gaus**2 * samp_freq/samples
'with a peak at 0 Hz with amplitude {} ({})').format(
expected[0], power[0]))
- if TEST_PLOTS :
+ if TEST_PLOTS:
figure = _pyplot.figure()
time_axes = figure.add_subplot(2, 1, 1)
time_axes.plot(
freq_axes.plot(freq_axis, expected, 'b-')
freq_axes.set_title('freq series')
-def _test_unitary_power_spectrum_gaussian_suite() :
+def _test_unitary_power_spectrum_gaussian_suite():
print('Test unitary power spectrums on various gaussian functions')
- _test_unitary_power_spectrum_gaussian(area=1, std=1, samp_freq=10.0, samples=1024)
- _test_unitary_power_spectrum_gaussian(area=1, std=2, samp_freq=10.0, samples=1024)
- _test_unitary_power_spectrum_gaussian(area=1, std=1, samp_freq=10.0, samples=2048)
- _test_unitary_power_spectrum_gaussian(area=1, std=1, samp_freq=20.0, samples=2048)
- _test_unitary_power_spectrum_gaussian(area=3, std=1, samp_freq=10.0, samples=1024)
+ _test_unitary_power_spectrum_gaussian(
+ area=1, std=1, samp_freq=10.0, samples=1024)
+ _test_unitary_power_spectrum_gaussian(
+ area=1, std=2, samp_freq=10.0, samples=1024)
+ _test_unitary_power_spectrum_gaussian(
+ area=1, std=1, samp_freq=10.0, samples=2048)
+ _test_unitary_power_spectrum_gaussian(
+ area=1, std=1, samp_freq=20.0, samples=2048)
+ _test_unitary_power_spectrum_gaussian(
+ area=3, std=1, samp_freq=10.0, samples=1024)
_test_unitary_power_spectrum_gaussian(
area=_numpy.pi, std=_numpy.sqrt(2), samp_freq=_numpy.exp(1),
samples=1024)
-def window_hann(length) :
+def window_hann(length):
"Returns a Hann window array with length entries"
win = _numpy.zeros((length,), dtype=_numpy.float)
- for i in range(length) :
+ for i in range(length):
win[i] = 0.5*(1.0-_numpy.cos(2.0*_numpy.pi*_numpy.float(i)/(length)))
# avg value of cos over a period is 0
# so average height of Hann window is 0.5
return win
def avg_power_spectrum(data, freq=1.0, chunk_size=2048,
- overlap=True, window=window_hann) :
+ overlap=True, window=window_hann):
"""Compute the avgerage power spectrum of DATA.
Taken with a sampling frequency FREQ.
'chunk_size {} should be a power of 2'.format(chunk_size))
nchunks = len(data)/chunk_size # integer division = implicit floor
- if overlap :
+ if overlap:
chunk_step = chunk_size/2
- else :
+ else:
chunk_step = chunk_size
win = window(chunk_size) # generate a window of the appropriate size
# >>> help(numpy.fft.fftpack.rfft) for Numpy's explaination.
# See Numerical Recipies for a details.
power = _numpy.zeros((chunk_size/2+1,), dtype=_numpy.float)
- for i in range(nchunks) :
+ for i in range(nchunks):
starti = i*chunk_step
stopi = starti+chunk_size
fft_chunk = _numpy.fft.rfft(data[starti:stopi]*win)
return (freq_axis, power)
def unitary_avg_power_spectrum(data, freq=1.0, chunk_size=2048,
- overlap=True, window=window_hann) :
+ overlap=True, window=window_hann):
"""Compute the average power spectrum, preserving normalization
"""
freq_axis,power = avg_power_spectrum(
data, freq, chunk_size, overlap, window)
- # 2.0 / (freq * chunk_size) |rfft()|**2 --> unitary_power_spectrum
+ # 2.0 / (freq * chunk_size) |rfft()|**2 --> unitary_power_spectrum
# see unitary_power_spectrum()
power *= 2.0 / (freq*_numpy.float(chunk_size)) * 8/3
- # * 8/3 to remove power from windowing
+ # * 8/3 to remove power from windowing
# <[x(t)*w(t)]**2> = <x(t)**2 * w(t)**2> ~= <x(t)**2> * <w(t)**2>
# where the ~= is because the frequency of x(t) >> the frequency of w(t).
# So our calulated power has and extra <w(t)**2> in it.
- # For the Hann window, <w(t)**2> = <0.5(1 + 2cos + cos**2)> = 1/4 + 0 + 1/8 = 3/8
- # For low frequency components, where the frequency of x(t) is ~= the frequency of w(t),
- # The normalization is not perfect. ??
+ # For the Hann window,
+ # <w(t)**2> = <0.5(1 + 2cos + cos**2)> = 1/4 + 0 + 1/8 = 3/8
+ # For low frequency components,
+ # where the frequency of x(t) is ~= the frequency of w(t),
+ # the normalization is not perfect. ??
# The normalization approaches perfection as chunk_size -> infinity.
return (freq_axis, power)
-def _test_unitary_avg_power_spectrum_sin(sin_freq=10, samp_freq=512, samples=1024,
- chunk_size=512, overlap=True,
- window=window_hann) :
+def _test_unitary_avg_power_spectrum_sin(
+ sin_freq=10, samp_freq=512, samples=1024, chunk_size=512, overlap=True,
+ window=window_hann):
x = _numpy.zeros((samples,), dtype=_numpy.float)
samp_freq = _numpy.float(samp_freq)
- for i in range(samples) :
+ for i in range(samples):
x[i] = _numpy.sin(2.0 * _numpy.pi * (i/samp_freq) * sin_freq)
freq_axis, power = unitary_avg_power_spectrum(
x, samp_freq, chunk_size, overlap, window)
sin_freq, expected[i], freq_axis[imax], power[imax]))
Pexp = 0
P = 0
- for i in range(len(freq_axis)) :
+ for i in range(len(freq_axis)):
Pexp += expected[i] * df
P += power[i] * df
print('The total power should be {} ({})'.format(Pexp, P))
- if TEST_PLOTS :
+ if TEST_PLOTS:
figure = _pyplot.figure()
time_axes = figure.add_subplot(2, 1, 1)
time_axes.plot(
freq_axes.set_title(
'{} samples of sin at {} Hz'.format(samples, sin_freq))
-def _test_unitary_avg_power_spectrum_sin_suite() :
+def _test_unitary_avg_power_spectrum_sin_suite():
print('Test unitary avg power spectrums on variously shaped sin functions')
- _test_unitary_avg_power_spectrum_sin(sin_freq=5, samp_freq=512, samples=1024)
- _test_unitary_avg_power_spectrum_sin(sin_freq=5, samp_freq=512, samples=2048)
- _test_unitary_avg_power_spectrum_sin(sin_freq=5, samp_freq=512, samples=4098)
- _test_unitary_avg_power_spectrum_sin(sin_freq=17, samp_freq=512, samples=1024)
- _test_unitary_avg_power_spectrum_sin(sin_freq=5, samp_freq=1024, samples=2048)
+ _test_unitary_avg_power_spectrum_sin(
+ sin_freq=5, samp_freq=512, samples=1024)
+ _test_unitary_avg_power_spectrum_sin(
+ sin_freq=5, samp_freq=512, samples=2048)
+ _test_unitary_avg_power_spectrum_sin(
+ sin_freq=5, samp_freq=512, samples=4098)
+ _test_unitary_avg_power_spectrum_sin(
+ sin_freq=17, samp_freq=512, samples=1024)
+ _test_unitary_avg_power_spectrum_sin(
+ sin_freq=5, samp_freq=1024, samples=2048)
# test long wavelenth sin, so be closer to window frequency
- _test_unitary_avg_power_spectrum_sin(sin_freq=1, samp_freq=1024, samples=2048)
- # finally, with some irrational numbers, to check that I'm not getting lucky
+ _test_unitary_avg_power_spectrum_sin(
+ sin_freq=1, samp_freq=1024, samples=2048)
+ # finally, with some irrational numbers, to check that I'm not
+ # getting lucky
_test_unitary_avg_power_spectrum_sin(
sin_freq=_numpy.pi, samp_freq=100*_numpy.exp(1), samples=1024)
-def test() :
+def test():
_test_rfft_suite()
_test_unitary_rfft_parsevals_suite()
_test_unitary_rfft_rect_suite()
projects / FFT-tools.git / blobdiff