From e746d9acaab757afbe480cc0d8befc71bca3ecf7 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Sun, 18 Nov 2012 17:19:42 -0500 Subject: [PATCH] FFT_tools: wrap long comments for PEP8 compliance --- FFT_tools.py | 39 +++++++++++++++++++++++++-------------- 1 file changed, 25 insertions(+), 14 deletions(-) diff --git a/FFT_tools.py b/FFT_tools.py index 352083a..87dcce5 100644 --- a/FFT_tools.py +++ b/FFT_tools.py @@ -296,17 +296,20 @@ def 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 @@ -386,7 +389,7 @@ def _test_unitary_power_spectrum_sin_suite(): _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 @@ -430,8 +433,10 @@ 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) @@ -453,9 +458,12 @@ def _test_unitary_power_spectrum_gaussian( # 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 @@ -553,16 +561,18 @@ def unitary_avg_power_spectrum(data, freq=1.0, chunk_size=2048, """ 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> = ~= * # where the ~= is because the frequency of x(t) >> the frequency of w(t). # So our calulated power has and extra in it. - # For the Hann window, = <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, + # = <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) @@ -612,7 +622,8 @@ def _test_unitary_avg_power_spectrum_sin_suite(): _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 + # 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) -- 2.26.2