Add upsampling via phase-tranformation experiment.

diff --git a/misc/resampling.py b/misc/resampling.py
index e3eb49c9ce623eaed67420b8f396d22a6e2298ca..59cd038811e8ed10b2793395c7e0d8662bc06364 100755 (executable)
--- a/misc/resampling.py
+++ b/misc/resampling.py
@@ -5,12 +5,38 @@
 # causes the non-upsampled integral to become invalid when the number of
 # samples is even.
 
-from numpy import concatenate, append, array, size, linalg
+from numpy import concatenate, append, array, size, linalg, pi, cos, sin, zeros, complex64
 from scipy import fft, ifft
 
 def upsample(signal, multiple):
   return ifft(pad(fft(signal), multiple))
 
+def upsample_phase_shift(signal):
+  transformed = fft(signal)
+  theta_base = pi/size(transformed)
+
+  for freq in range(size(transformed)):
+    theta = -theta_base * freq
+
+    if size(transformed) % 2 == 0 and freq == (size(transformed) / 2):
+      phi = 0
+    else:
+      if (freq < 1 + size(transformed)/2):
+        theta = theta_base * freq
+      else:
+        theta = pi + (theta_base * freq)
+      phi = cos(theta) + sin(theta)*1j
+
+    transformed[freq] *= phi
+
+  shifted_signal = ifft(transformed)
+
+  output = zeros(size(signal)*2, complex64)
+  output[0::2] = signal
+  output[1::2] = shifted_signal
+
+  return output
+
 def pad(signal, multiple):
   split = size(signal) / 2
   if (size(signal) % 2 == 0):
@@ -51,6 +77,14 @@ def compareIntegrals(signal):
   print "Upsampled 4x integral: " + str(integralVeryFine.round(5))
   print "Absolute difference: " + str(abs(integralCoarse-integralVeryFine))
 
+def compareUpsampling(signal):
+  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))
+
+
 def checkUpsampling(signal):
   upsampled = upsample(signal, 4)
   decimated = everyN(upsampled, 4)
@@ -67,6 +101,10 @@ checkUpsampling(signal)
 print
 checkUpsampling(signal2)
 print "\n"
+compareUpsampling(signal)
+print "\n"
+compareUpsampling(signal2)
+print "\n"
 compareIntegrals(signal)
 print "\n"
 compareIntegrals(signal2)