From b031c3d1bee4b0405d2b5f92606fa5760559196d Mon Sep 17 00:00:00 2001 From: Francis Russell Date: Tue, 26 Mar 2013 17:22:27 +0000 Subject: [PATCH] Reformat resampling.py. --- misc/resampling.py | 137 +++++++++++++++++++++++++-------------------- 1 file changed, 75 insertions(+), 62 deletions(-) diff --git a/misc/resampling.py b/misc/resampling.py index 59cd038..7c74320 100755 --- a/misc/resampling.py +++ b/misc/resampling.py @@ -5,93 +5,104 @@ # 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))) -- 2.47.3