1 #!/usr/bin/python 2 3 # Copyright (C) 2012 The Android Open Source Project 4 # 5 # Licensed under the Apache License, Version 2.0 (the "License"); 6 # you may not use this file except in compliance with the License. 7 # You may obtain a copy of the License at 8 # 9 # http://www.apache.org/licenses/LICENSE-2.0 10 # 11 # Unless required by applicable law or agreed to in writing, software 12 # distributed under the License is distributed on an "AS IS" BASIS, 13 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 # See the License for the specific language governing permissions and 15 # limitations under the License. 16 17 import numpy as np 18 import numpy.linalg 19 import scipy as sp 20 import scipy.fftpack 21 import scipy.signal 22 import math 23 import sys 24 from multiprocessing import Pool 25 26 def convolution(data0, data1reversed, n): 27 """calculate convolution part of data0 with data1 from pos n""" 28 N = len(data1reversed) 29 return np.dot(data0[n:N+n], data1reversed) 30 31 32 def convolutionstar(args): 33 return convolution(*args) 34 35 def calc_delay(data0, data1): 36 """Calcuate delay between two data. data0 is assumed to be recorded first, 37 and will have longer length than data1 38 returns delay between data0 and data1 in number of samples in data0's point of view""" 39 40 len0 = len(data0) 41 len1 = len(data1) 42 if len1 > len0: 43 print "data1 longer than data0" 44 return -1 45 searchLen = len0 - len1 46 data1reverse = data1[::-1] 47 48 # This is faster than signal.correlate as there is no need to process 49 # full data, but still it is slow. about 18 secs for data0 of 4 secs with data1 of 1 secs 50 print "***Caluclating delay, may take some time***" 51 gData0 = data0 52 gData1 = data1reverse 53 pool = Pool(processes = 4) 54 TASK = [(data0, data1reverse, i) for i in range(searchLen)] 55 result = pool.map(convolutionstar, TASK) 56 57 return np.argmax(result) 58 59 60 # test code 61 if __name__=="__main__": 62 samplingRate = 44100 63 durationInSec = 0.001 64 if len(sys.argv) > 1: 65 durationInSec = float(sys.argv[1]) 66 signalFrequency = 1000 67 samples = float(samplingRate) * float(durationInSec) 68 index = np.linspace(0.0, samples, num=samples, endpoint=False) 69 time = index / samplingRate 70 multiplier = 2.0 * np.pi * signalFrequency / float(samplingRate) 71 data0 = np.sin(index * multiplier) 72 DELAY = durationInSec / 2.0 * samplingRate 73 data1 = data0[DELAY:] 74 delay = calc_delay(data0, data1) 75 print "calc_delay returned", delay, " while expecting ", DELAY 76 77 78
