Add experiment on validity of non-upsampled overlap integrals.

diff --git a/.gitignore b/.gitignore
index 7d33312aa21d6b225ec2b1a99f8772ab24fe5b1c..4ef2c677171de01ac4a7396368e9331bcd03acb2 100644 (file)
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 /dist/
 ofc.jar
 .*.swp
+*.pyc

diff --git a/misc/resampling.py b/misc/resampling.py
new file mode 100755 (executable)
index 0000000..40e7701
--- /dev/null
+++ b/misc/resampling.py
@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+
+# Evaluates the integral of a squared signal with and without resampling
+# to a higher number of samples. The result is suggestive that aliasing
+# causes the non-upsampled integral to become invalid when the number of
+# samples is even.
+
+from numpy import concatenate, append, array, size, linalg
+from scipy import fft, ifft
+
+def upsample(signal, multiple):
+  return ifft(pad(fft(signal), multiple))
+
+def pad(signal, multiple):
+  split = size(signal) / 2
+  if (size(signal) % 2 == 0):
+    nyquist = signal[split]
+    upsampled = concatenate(
+      [
+        signal[:split],
+        [nyquist/2.0],
+        [0]*(size(signal)*(multiple-1)-1),
+        [nyquist/2.0],
+        signal[split+1:]
+      ])
+  else:
+    upsampled = concatenate(
+      [
+        signal[:split+1],
+        [0]*size(signal)*(multiple-1),
+        signal[split+1:]
+      ])
+
+  return upsampled*multiple
+
+def everyN(signal, n):
+  return array([signal[n*i] for i in range((1 + size(signal))/n)])
+
+def integrate(x):
+  return sum(x)/float(size(x))
+
+def compareIntegrals(signal):
+  integralCoarse = integrate(signal**2)
+  integralFine = integrate(upsample(signal, 2)**2)
+  integralVeryFine = integrate(upsample(signal, 4)**2)
+
+  print "Signal: " + str(signal)
+  print "Signal length: " + str(size(signal)) + " samples"
+  print "Integral: " + str(integralCoarse.round(5))
+  print "Upsampled 2x integral: " + str(integralFine.round(5))
+  print "Upsampled 4x integral: " + str(integralVeryFine.round(5))
+  print "Absolute difference: " + str(abs(integralCoarse-integralVeryFine))
+
+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))
+
+signal = array([1,3.5,-7,2.8,51,-3,-23,1000])
+signal2 = append(signal, [1])
+
+checkUpsampling(signal)
+print
+checkUpsampling(signal2)
+print "\n"
+compareIntegrals(signal)
+print "\n"
+compareIntegrals(signal2)
+