From 144fd1178e0275237b80879e8fb6859296a4c1fa Mon Sep 17 00:00:00 2001
From: Francis Russell <francis@unchartedbackwaters.co.uk>
Date: Tue, 9 Apr 2013 22:23:50 +0100
Subject: [PATCH] Compare upsampling against generated reference signal.

---
 misc/resampling.py | 64 ++++++++++++++++++++++++++++------------------
 1 file changed, 39 insertions(+), 25 deletions(-)

diff --git a/misc/resampling.py b/misc/resampling.py
index ff6fa8c..de96488 100755
--- a/misc/resampling.py
+++ b/misc/resampling.py
@@ -6,7 +6,7 @@
 # samples is even.
 
 from numpy import (concatenate, append, array, size, linalg, pi, cos, sin,
-    zeros, complex128)
+    zeros, complex128, linspace)
 from scipy import fft, ifft
 
 
@@ -14,6 +14,15 @@ def upsample(signal, multiple):
     return ifft(pad(fft(signal), multiple))
 
 
+def gensignal(components, num_samples):
+    space = linspace(0, 1, num_samples + 1)[0:num_samples]
+    signal = space * 0
+    for (freq, phase, magnitude) in components:
+        assert(freq <= num_samples / 2)
+        signal = signal + magnitude * sin(2 * pi * space * freq + phase)
+    return signal
+
+
 def upsample_phase_shift(signal):
     transformed = fft(signal)
     theta_base = pi / size(transformed)
@@ -65,7 +74,7 @@ def pad(signal, multiple):
 
 
 def everyN(signal, n):
-    return array([signal[n * i] for i in range((1 + size(signal)) / n)])
+    return signal[::n]
 
 
 def integrate(x):
@@ -85,40 +94,45 @@ def compareIntegrals(signal):
     print "Absolute difference: " + str(abs(integralCoarse - integralVeryFine))
 
 
-def compareUpsampling(signal):
+def compareUpsampling(signal_fine):
+    signal = everyN(signal_fine, 2)
     upsampled_padded = upsample(signal, 2)
     upsampled_rotated = upsample_phase_shift(signal)
 
     print ("Comparing upsampling using padding and rotation of signal: " +
         str(signal))
-    print ("L2-norm between corresponding samples: " +
-        str(linalg.norm(upsampled_padded - upsampled_rotated)))
-
+    print ("L2-norm between original and interpolated (padded): " +
+        str(linalg.norm(signal_fine - upsampled_padded)))
+    print ("L2-norm between original and interpolated (rotated): " +
+        str(linalg.norm(signal_fine - upsampled_rotated)))
+
+
+bandlimited1_fine = gensignal(
+        [(0, pi / 2, 10),
+         (3, pi / 5, 7),
+         (5, pi / 2, 17)
+        ], 20)
+bandlimited2_fine = gensignal(
+        [(0, pi / 2, 7),
+         (3, pi / 5, 13),
+         (4, pi / 2, 11)
+        ], 18)
+
+compareUpsampling(bandlimited1_fine)
+print "\n"
+compareUpsampling(bandlimited2_fine)
+print "\n"
 
-def checkUpsampling(signal):
-    upsampled = upsample(signal, 4)
-    decimated = everyN(upsampled, 4)
-    print "Comparing original and upsampled version of signal: " + str(signal)
-    print ("L2-norm between corresponding samples: " +
-        str(linalg.norm(decimated - signal)))
+random1 = array([1, 3.5, -7, 2.8, 51, -3, -23, 1000])
+random2 = append(random1, [1])
 
-signal = array([1, 3.5, -7, 2.8, 51, -3, -23, 1000])
-signal2 = append(signal, [1])
+compareIntegrals(random1)
+print "\n"
+compareIntegrals(random2)
 
 evenSignalFreqs = array([0., 0., 0., 0., 1., 0., 0., 0.])
 oddSignalFreqs = array([0., 0., 0., 0.5, 0.5, 0., 0.])
 
-checkUpsampling(signal)
-print
-checkUpsampling(signal2)
-print "\n"
-compareUpsampling(signal)
-print "\n"
-compareUpsampling(signal2)
-print "\n"
-compareIntegrals(signal)
-print "\n"
-compareIntegrals(signal2)
 print "\nSquaring signal with even number of samples..."
 print "Original signal (in frequency domain): " + str(evenSignalFreqs)
 print ("Squared signal (in frequency domain): " +
-- 
2.47.3