]> git.unchartedbackwaters.co.uk Git - francis/ofc.git/commitdiff
Reformat resampling.py.
authorFrancis Russell <francis@unchartedbackwaters.co.uk>
Tue, 26 Mar 2013 17:22:27 +0000 (17:22 +0000)
committerFrancis Russell <francis@unchartedbackwaters.co.uk>
Tue, 26 Mar 2013 17:37:19 +0000 (17:37 +0000)
misc/resampling.py

index 59cd038811e8ed10b2793395c7e0d8662bc06364..7c74320f3acbca86a028e1b55a4b4d8830fe7c96 100755 (executable)
 # causes the non-upsampled integral to become invalid when the number of
 # samples is even.
 
-from numpy import concatenate, append, array, size, linalg, pi, cos, sin, zeros, complex64
+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))
+    return ifft(pad(fft(signal), multiple))
+
 
 def upsample_phase_shift(signal):
-  transformed = fft(signal)
-  theta_base = pi/size(transformed)
+    transformed = fft(signal)
+    theta_base = pi / size(transformed)
 
-  for freq in range(size(transformed)):
-    theta = -theta_base * freq
+    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
+        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
+        transformed[freq] *= phi
 
-  shifted_signal = ifft(transformed)
+    shifted_signal = ifft(transformed)
 
-  output = zeros(size(signal)*2, complex64)
-  output[0::2] = signal
-  output[1::2] = shifted_signal
+    output = zeros(size(signal) * 2, complex64)
+    output[0::2] = signal
+    output[1::2] = shifted_signal
+
+    return output
 
-  return output
 
 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
+    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)])
+    return array([signal[n * i] for i in range((1 + size(signal)) / n)])
+
 
 def integrate(x):
-  return sum(x)/float(size(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)
+    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))
 
-  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 compareUpsampling(signal):
-  upsampled_padded = upsample(signal, 2)
-  upsampled_rotated = upsample_phase_shift(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))
+    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)
-  print "Comparing original and upsampled version of signal: " + str(signal)
-  print "L2-norm between corresponding samples: " + str(linalg.norm(decimated - 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])
+signal = array([1, 3.5, -7, 2.8, 51, -3, -23, 1000])
 signal2 = append(signal, [1])
 
 evenSignalFreqs = array([0., 0., 0., 0., 1., 0., 0., 0.])
@@ -110,8 +121,10 @@ 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): " + str(fft(size(evenSignalFreqs)*ifft(evenSignalFreqs)**2).round(5))
+print ("Squared signal (in frequency domain): " +
+    str(fft(size(evenSignalFreqs) * ifft(evenSignalFreqs) ** 2).round(5)))
 
 print "\nSquaring signal with odd number of samples..."
 print "Original signal (in frequency domain): " + str(oddSignalFreqs)
-print "Squared signal (in frequency domain): " + str(fft(size(oddSignalFreqs)*ifft(oddSignalFreqs)**2).round(5))
+print ("Squared signal (in frequency domain): " +
+    str(fft(size(oddSignalFreqs) * ifft(oddSignalFreqs) ** 2).round(5)))