From 2229b7919df93e1b29af4e7617e95a195acf64c8 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Sat, 4 Oct 2008 09:43:54 -0400 Subject: [PATCH] Got Parseval's checks working. --- FFT_tools.py | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/FFT_tools.py b/FFT_tools.py index 4abb0a7..b6b7b35 100644 --- a/FFT_tools.py +++ b/FFT_tools.py @@ -54,17 +54,18 @@ def _test_rfft(xs, Xs) : j = complex(0,1) n = len(xs) Xa = [] - for k in range(len(Xs)) : + for k in range(n) : Xa.append(sum([x*exp(-j*2*pi*k*m/n) for x,m in zip(xs,range(n))])) - assert (Xs[k]-Xa[k])/abs(Xa[k]) < 1e-6, \ - "rfft mismatch on element %d: %g != %g, relative error %g" \ - % (k, Xs[k], Xa[k], (Xs[k]-Xa[k])/abs(Xa[k])) + if k < len(Xs): + assert (Xs[k]-Xa[k])/abs(Xa[k]) < 1e-6, \ + "rfft mismatch on element %d: %g != %g, relative error %g" \ + % (k, Xs[k], Xa[k], (Xs[k]-Xa[k])/abs(Xa[k])) # Which should satisfy the discrete form of Parseval's theorem # n-1 n-1 # SUM |x_m|^2 = 1/n SUM |X_k|^2. # m=0 k=0 timeSum = sum([abs(x)**2 for x in xs]) - freqSum = sum([abs(X)**2 for X in Xa]) + freqSum = sum([abs(X)**2 for X in Xa]) assert abs(freqSum/float(n) - timeSum)/timeSum < 1e-6, \ "Mismatch on Parseval's, %g != 1/%d * %g" % (timeSum, n, freqSum) @@ -118,23 +119,31 @@ def unitary_rfft(data, freq=1.0) : freq_axis = linspace(0, freq/2, nsamps/2+1) return (freq_axis, trans) -def _test_unitary_rfft_parsevals(): +def _test_unitary_rfft_parsevals(xs, freq, freqs, Xs): 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] - dt = pi - freqs,Xs = unitary_rfft(xs, 1.0/dt) # Which should satisfy the discretized integral form of Parseval's theorem # n-1 n-1 # SUM |x_m|^2 dt = SUM |X_k|^2 df # m=0 k=0 + dt = 1.0/freq df = freqs[1]-freqs[0] assert (df - 1/(len(xs)*dt))/df < 1e-6, \ "Mismatch in spacing, %g != 1/(%d*%g)" % (df, len(xs), dt) + Xa = list(Xs) + for k in range(len(Xs)-1,1,-1) : + Xa.append(Xa[k]) + assert len(xs) == len(Xa), "Length mismatch %d != %d" % (len(xs), len(Xa)) lhs = sum([abs(x)**2 for x in xs]) * dt - rhs = sum([abs(X)**2 for X in Xs]) * df - assert abs(lhs - rhs)/lhs < 1e-6, "Mismatch on Parseval's, %g != %g" \ + rhs = sum([abs(X)**2 for X in Xa]) * df + assert abs(lhs - rhs)/lhs < 1e-4, "Mismatch on Parseval's, %g != %g" \ % (lhs, rhs) +def _test_unitary_rfft_parsevals_suite(): + xs = [1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1] + dt = pi + freqs,Xs = unitary_rfft(xs, 1.0/dt) + _test_unitary_rfft_parsevals(xs, 1.0/dt, freqs, Xs) + def _rect(t) : if abs(t) < 0.5 : return 1 @@ -556,7 +565,7 @@ def _test_unitary_avg_power_spectrum_sin_suite() : def test() : _test_rfft_suite() - _test_unitary_rfft_parsevals() + _test_unitary_rfft_parsevals_suite() _test_unitary_rfft_rect_suite() _test_unitary_rfft_gaussian_suite() _test_unitary_power_spectrum_sin_suite() -- 2.26.2