Got Parseval's checks working.

diff --git a/FFT_tools.py b/FFT_tools.py
index 4abb0a7523ad47b3479c0910fd41a4f35d0fe2d5..b6b7b353ecf12dd25c9ed1da8be136ce50d223d8 100644 (file)
--- 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()