From: Francis Russell <francis@unchartedbackwaters.co.uk>
Date: Tue, 26 Mar 2013 15:57:46 +0000 (+0000)
Subject: Add upsampling via phase-tranformation experiment.
X-Git-Url: https://git.unchartedbackwaters.co.uk/w/?a=commitdiff_plain;h=aff04a7db8d0cf4b94979f864b30c798364b34bd;p=francis%2Fofc.git

Add upsampling via phase-tranformation experiment.
---

diff --git a/misc/resampling.py b/misc/resampling.py
index e3eb49c..59cd038 100755
--- 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)