num = 2**_numpy.floor(lnum)
return 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)
num = 2**_numpy.round(lnum)
return num
+
def ceil_pow_of_two(num):
"Round num up to the closest exact a power of two."
lnum = _numpy.log2(num)
num = 2**_numpy.ceil(lnum)
return num
+
def _test_rfft(xs, Xs):
# Numpy's FFT algoritm returns
# n-1
"Mismatch on Parseval's, {} != 1/{} * {}".format(
timeSum, n, freqSum))
+
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):
"""Compute the Fourier transform of real data.
freq_axis = _numpy.linspace(0, freq/2, nsamps/2+1)
return (freq_axis, trans)
+
def _test_unitary_rfft_parsevals(xs, freq, freqs, Xs):
# Which should satisfy the discretized integral form of Parseval's theorem
# n-1 n-1
if _numpy.abs(lhs - rhs)/lhs >= 1e-4:
raise ValueError("Mismatch on Parseval's, {} != {}".format(lhs, rhs))
+
def _test_unitary_rfft_parsevals_suite():
print("Test unitary rfft on Parseval's theorem")
xs = [1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1]
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:
return 1
else:
return 0
+
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)"
freq_axes.plot(freq_axis, expected, 'b-')
freq_axes.set_title('freq series')
+
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=0.7, samp_freq=50, samples=512)
_test_unitary_rfft_rect(a=3.0, samp_freq=60, samples=1024)
+
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):
"Show fft(rect(at)) = 1/abs(a) * sinc(f/a)"
freq_axes.plot(freq_axis, expected, 'b-')
freq_axes.set_title('freq series')
+
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=3.0, samp_freq=60, samples=1024)
-
def power_spectrum(data, freq=1.0):
"""Compute the power spectrum of DATA taken with a sampling frequency FREQ.
power = (trans * trans.conj()).real # We want the square of the amplitude.
return (freq_axis, power)
+
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
return (freq_axis, power)
+
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)
freq_axes.set_title(
'{} samples of sin at {} Hz'.format(samples, sin_freq))
+
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 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):
x = _numpy.zeros((samples,), dtype=_numpy.float)
samp_freq = _numpy.float(samp_freq)
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():
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=_numpy.pi, samp_freq=_numpy.exp(1), samples=1024)
+
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):
x = _numpy.zeros((samples,), dtype=_numpy.float)
freq_axes.plot(freq_axis, expected, 'b-')
freq_axes.set_title('freq series')
+
def _test_unitary_power_spectrum_gaussian_suite():
print('Test unitary power spectrums on various gaussian functions')
_test_unitary_power_spectrum_gaussian(
area=_numpy.pi, std=_numpy.sqrt(2), samp_freq=_numpy.exp(1),
samples=1024)
+
def window_hann(length):
"Returns a Hann window array with length entries"
win = _numpy.zeros((length,), dtype=_numpy.float)
# 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):
"""Compute the avgerage power spectrum of DATA.
power /= _numpy.float(nchunks)
return (freq_axis, power)
+
def unitary_avg_power_spectrum(data, freq=1.0, chunk_size=2048,
overlap=True, window=window_hann):
"""Compute the average power spectrum, preserving normalization
# 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):
freq_axes.set_title(
'{} samples of sin at {} Hz'.format(samples, sin_freq))
+
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(
_test_unitary_power_spectrum_gaussian_suite()
_test_unitary_avg_power_spectrum_sin_suite()
+
if __name__ == '__main__':
from optparse import OptionParser
projects / FFT-tools.git / blobdiff