# 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.])
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)))