From: W. Trevor King <wking@drexel.edu>
Date: Sat, 4 Oct 2008 13:43:54 +0000 (-0400)
Subject: Got Parseval's checks working.
X-Git-Tag: v0.2~5
X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=2229b7919df93e1b29af4e7617e95a195acf64c8;p=stepper.git

Got Parseval's checks working.
---

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()