Source Code for Module wiat.Data

  1  # Copyright 2008 Antoine Lefebvre 
  2  # 
  3  # This file is part of "The Wind Instrument Acoustic Toolkit". 
  4  # 
  5  # "The Wind Instrument Acoustic Toolkit" is free software: 
  6  # you can redistribute it and/or modify it under the terms of 
  7  # the GNU General Public License as published by the 
  8  # Free Software Foundation, either version 3 of the License, or 
  9  # (at your option) any later version. 
 10  # 
 11  # "The Wind Instrument Acoustic Toolkit" is distributed in the hope 
 12  # that it will be useful, but WITHOUT ANY WARRANTY; without even 
 13  # the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
 14  # See the GNU General Public License for more details. 
 15  # 
 16  # You should have received a copy of the GNU General Public License 
 17  # along with "The Wind Instrument Acoustic Toolkit". 
 18  # If not, see <http://www.gnu.org/licenses/>. 
 19  # 
 20  # All rights reserved. 
 21   
 22  """ Data.py - Functions to load data from Matlab files and to  
 23      interpolate the frequency and magnitude of the maxima. 
 24   
 25  """ 
 26   
 27  __all__ = ['LoadMatlabData', 'InterpolateMaxima'] 
 28   
 29  import scipy.io 
 30  from scipy import optimize, linalg 
 31  import numpy 
 32  from numpy import sign, diff 
 33   
 34  from wiat.util import calculate_interval_in_cents 
 35   
 36 -def load_from_mat(names, filename): 
 37      """ Load data from a matlab file. All variable name in the names tuple are  
 38          fetched and returned . 
 39   
 40          Example:: 
 41   
 42          >>> import os 
 43          >>> import numpy 
 44          >>> import wiat 
 45          >>> f,z = wiat.Data.load_from_mat(('freqs','z'), 
 46          ...    os.path.join(wiat.MEASUREMENTDIR,'20080412Ztoneholeless.mat')) 
 47          >>> numpy.size(f) == numpy.size(z) 
 48          True 
 49           
 50   
 51          @param names: A tuple containing variable names to get from the Matlab 
 52              file. 
 53          @param filename: the path to the Matlab file 
 54          @return: a list of the extracted arrays. None if the name do not exist. 
 55      """ 
 56   
 57      try: 
 58          data = scipy.io.loadmat(filename) 
 59      except IOError, arg: 
 60          print arg 
 61          raise IOError(arg) 
 62   
 63      ret = [] 
 64      for name in names: 
 65          if not data.has_key(name): 
 66              print 'Warning: no such variable %s in Matlab file %s'\ 
 67                      %(name,filename) 
 68          ret.append(data.get(name)) 
 69      return ret 
 70   
 71 -class __Parabola: 
 72      " Class for fiting an equation of the form M{y = ax**2 + bx + c} " 
 73 -    def __init__(self): 
 74          self.n = 0 
 75           
 76 -    def __call__(self, f): 
 77          c = self.C 
 78          return c[0]*f**2+c[1]*f+c[2] 
 79       
 80 -    def fit_points(self, px, py): 
 81          """ Fit the parabola passing through 3 points with 3 equations 
 82   
 83              M{y1 = a * x1**2 + b*x1 + c} 
 84              M{y2 = a * x2**2 + b*x2 + c} 
 85              M{y3 = a * x3**2 + b*x3 + c} 
 86   
 87          """ 
 88          self.px = px 
 89          self.py = py 
 90          B = py 
 91          tmp = numpy.array([px**2, px, [1.,1.,1.]]) 
 92          A = numpy.array([tmp[:,0], tmp[:,1], tmp[:,2]]) 
 93          self.C = scipy.linalg.solve(A, B) 
 94   
 95 -    def lowf(self): 
 96          return self.px[0] 
 97   
 98 -    def highf(self): 
 99          return self.px[2] 
100   
101 -    def xc(self): 
102          """ position of vertical axis of symmetry: h = -b / 2a """ 
103          return - self.C[1] / (2*self.C[0]) 
104   
105 -    def yc(self): 
106          """ value of parabola on axis of symmetry: k = (4ac - b**2) / 4a """ 
107          return (4*self.C[0]*self.C[2] - self.C[1]**2) / (4 * self.C[0]) 
108       
109 -def interpolate_maxima(f, Z, number_of_maxima=None, fmin=100.): 
110      """ Find the maxima in the data and interpolate with Parbola to find the  
111          frequency and amplitude of the maxima. 
112   
113          Results that are noisy are likely to pose problems to the algorithm. 
114   
115          @param fmin: minimal frequency we are interested in, useful to avoid  
116              low frequency noise to infect the results. 
117      """ 
118   
119      derivative = diff(abs(Z)) 
120      signarray = sign(derivative) 
121      crossing  = diff(signarray) 
122      indices   = crossing.nonzero() 
123   
124      maximalist = [] 
125   
126      for i in indices[0]: 
127          if f[i] < fmin: 
128              continue 
129   
130         # when there is noise in the data, some maxima appears at various  
131         # points in the curve. 
132         # TODO: find a method to eliminate those false maxima 
133   
134              # in order to eliminate them, we can check various criteria. 
135              # generally, the derivative of the phase is large at the maxima 
136          #if abs(angle(Z[i])-angle(Z[i+1]))/(f[i+1]-f[1]) < 0.001: 
137          #    print 'skipping f = %f'%f[i] 
138          #    continue 
139   
140          # minima are much smaller than 1, maxima are much larger 
141          if abs(Z[i]) > 1.: 
142              # find the indice of the maxima around i 
143              while True: 
144                  if abs(Z[i-1]) > abs(Z[i]): 
145                      i = i - 1 
146                  elif abs(Z[i+1]) > abs(Z[i]): 
147                      i = i + 1 
148                  else: 
149                      break 
150                       
151              P = __Parabola() 
152              P.fit_points(f[i-1:i+2], abs(Z[i-1:i+2])) 
153   
154              # if the new maxima is closer than a semitone to the previous one 
155              # it is likely a duplicate due to noise... the highest one is kept 
156              if maximalist != []: 
157                  if abs(calculate_interval_in_cents(P.xc(), maximalist[-1][0]))\ 
158                     < 300.: 
159                      if P.yc() > maximalist[-1][1]: 
160                          maximalist.pop() 
161                      else: 
162                          continue 
163                           
164              #parabolas.append([i, P])                 
165              maximalist.append( ( P.xc(), P.yc() ) ) 
166      return maximalist 
167   
168   
169  # Code to be adapted 
170  # 
171  #def fitfunc(p,f): 
172  #    return  ( p[0] / (1 + p[1]*(f - p[2])**2 ) ) 
173  # 
174  #def errfunc(p,f,y): 
175  #    ff = fitfunc(p,f) 
176  #    return (ff - y) 
177  # 
178  #def fitLorentzPeak(self): 
179  #    l = self.listMaxima() 
180  #    print "Fitting peaks for %s" % self.desc.object 
181  # 
182  #    self.LorentzPeaks = [] 
183  #         
184  #    for k in arange(0,len(l),1): 
185  #        i = self.maximalist[k] 
186  #        # test if it is a maximum or minimum and proceed for maximum only 
187  #        #print abs(self.Z[i]) 
188  #        if abs(self.Z[i]) > 1.5: 
189  #            #print 'working on maxima no. %d'%k 
190  #            # p0 contains three value, the peak height, a value related to  
191               # peak width, and frequency 
192  #            A = abs(self.Z[i]) 
193  #            p0 = [A, 0.05, self.f[i]]  # Initial guess for the parameters 
194  #            deltaf = self.f[1]-self.f[0] 
195  #            ni = int(self.f[i]*(pow(2.,1./48.)-1.)/deltaf)  
196                  # use a small interval around the center 
197  #            if ni <= 1: 
198  #                ni = 2 
199  #                     
200  #            nc = self.maximalist[k] 
201  #            #print ni, nc 
202  #            range = arange(nc-ni,nc+ni+1,1) 
203  # 
204  #            #print 'optimizing with points',  
205                #(self.f[range], abs(self.Z[range])) 
206  #                                               
207  #            p1,success = scipy.optimize.leastsq(errfunc, p0[:], 
208  #args = (self.f[range], abs(self.Z[range]))) 
209  #            #print p0, p1, success 
210  # 
211  #            f1 = p1[2] 
212  #            nc = int(round((f1-self.f[0])/deltaf)) 
213  #            range = arange(nc-ni,nc+ni,1) 
214  #            p1,success = scipy.optimize.leastsq(errfunc, p0[:],  
215  #args = (self.f[range], abs(self.Z[range]))) 
216  #            #print p0, p1, success 
217  #            # TODO: recalculate the center indice from p1[2]  
218  #(the frequency) and reoptimize 
219  #                
220  #            fmax=self.f[nc+ni] 
221  #            fmin=self.f[nc-ni] 
222  #            if success and p1[2] > 0.: 
223  #                self.LorentzPeaks.append([p1,fmin,fmax]) 
224  #            else: 
225  #                print 'maxima rejected' 
226