Source Code for Module wiat.TM.Inst

  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  """ Inst.py - This file contains the class defining an Instrument and the 
 23      functions CalculateImpedance and CalculateMaxima. 
 24   
 25      The CalculateImpedance returns the impedance of an instrument for each 
 26      requested frequency. 
 27   
 28      The CalculateMaxima seach for the exact location of the maxima of the 
 29      instrument and return a list of (frequency,magnitude) tuples. 
 30           
 31  """ 
 32   
 33  __all__ = ['Instrument','CalculateImpedance','CalculateMaxima'] 
 34   
 35  import numpy 
 36  import scipy.optimize 
 37   
 38  from Waveguides import ToneHole 
 39   
 40 -def create_identity_matrix_array(n): 
 41      """ Returns an array of 2X2 identity matrices for uses in the transmission 
 42          matrix calculation method. 
 43      """ 
 44      M = numpy.array([[numpy.ones(n),numpy.zeros(n)], 
 45                       [numpy.zeros(n),numpy.ones(n)]]) 
 46      return M 
 47   
 48 -def multiply_matrix_arrays(TM1,TM2): 
 49      """ Multiplication of two arrays of 2X2 matrices. 
 50   
 51          Example:: 
 52           
 53          >>> TM1 = create_identity_matrix_array(10) 
 54          >>> TM2 = multiply_matrix_arrays(TM1,TM1) 
 55          >>> TM1 == TM2 
 56          array([[[ True,  True,  True,  True,  True,  True,  True,  True,  True, 
 57                    True], 
 58                  [ True,  True,  True,  True,  True,  True,  True,  True,  True, 
 59                    True]], 
 60          <BLANKLINE> 
 61                 [[ True,  True,  True,  True,  True,  True,  True,  True,  True, 
 62                    True], 
 63                  [ True,  True,  True,  True,  True,  True,  True,  True,  True, 
 64                    True]]], dtype=bool) 
 65          >>> 
 66   
 67      """ 
 68      n = range(numpy.shape(TM1)[-1]) 
 69      for i in n: 
 70          if i == 0: 
 71              R = numpy.dot(TM1[...,i],TM2[...,i]) 
 72          else: 
 73              T = numpy.dot(TM1[...,i],TM2[...,i]) 
 74              R = numpy.dstack((R, T)) 
 75      if numpy.size(n) == 1: 
 76          return numpy.resize(R,(2,numpy.shape(TM2)[1],1)) 
 77      else: 
 78          return R 
 79   
 80           
 81 -class Instrument: 
 82      """ Class used to descibed the geometry of a musical instrument. 
 83      """ 
 84 -    def __init__(self, name, sections, radiation, 
 85                   toneholes=None, fingerings=None): 
 86          """ Instrument class 
 87   
 88              Example:: 
 89   
 90              >>> import wiat 
 91              >>> from wiat.TM import * 
 92              >>> I1 = Instrument("Cylinder with one tone hole", 
 93              ...                [Cylinder(0.014,0.583)], UnflangedOpenEnd(0.014), 
 94              ...                [ToneHole(0.0117,0.0043,0.474-0.0117/2), 
 95              ...                 ToneHole(0.0112,0.0043,0.505-0.0112/2)], 
 96              ...                {'D':[1,1], 'C sharp':[0,1], 'C':[0, 0]}) 
 97   
 98              By default, all toneholes of the instrument are supposed close. Then 
 99              the fingering can be set with the set_fingering method. 
100   
101   
102              @param name: the name of the instrument 
103              @type name: string 
104              @param sections: a list of Waveguides 
105              @type sections: list 
106              @param radiation: a radiation condition 
107              @type radiation: an instance of a radiation class (see Radiation.py) 
108              @param toneholes: optional list of ToneHole  
109              @type toneholes: list 
110              @param fingerings: an optional fingering dictionary 
111                  The fingering dictionary should contains lists of open and  
112                  closed toneholes. The length of the lists should be the same  
113                  as the number of toneholes. 
114              @type fingerings: dictionary 
115   
116          """ 
117          self.name = name 
118          assert type(sections) is type([]), \ 
119              'The second argument should be a list.' 
120          self.sections = sections 
121          self.radiation = radiation 
122          self.toneholes = toneholes 
123   
124          # self.fingerings is a dictionary 
125          # self.fingering is a list 
126   
127          self.fingerings = fingerings 
128          if toneholes is not None: 
129              self.fingering = [0]*len(toneholes) # default closed toneholes 
130   
131 -    def __str__(self): 
132          string = '==============================================\n' 
133          string += "Instrument: %s\n\n"%self.name 
134   
135          string += "list of sections\n" 
136          string += "----------------\n" 
137          for s in self.sections: 
138              string = string + str(s) + '\n' 
139           
140          string += '\nradiation\n' 
141          string += '----------\n' 
142          string += str(self.radiation)+'\n' 
143          string += '==============================================\n' 
144          return string 
145   
146 -    def add_tonehole(self,tonehole,state=1): 
147          """ Add a new tonehole to the instrument. 
148   
149          This method is useful in the design stage of an instrument. 
150   
151          @param tonehole: the new tonehole to be added. 
152          @type tonehole: ToneHole 
153          @param state: if the tonehole is to be added in the open (1) or  
154          closed (0) state. 
155          """ 
156          if not isinstance(tonehole,ToneHole): 
157              raise TypeError('Expecting a tonehole') 
158   
159          if self.toneholes is None: 
160              self.toneholes =  [tonehole] 
161              self.fingering =  [state] 
162          else: 
163              self.toneholes.append(tonehole) 
164              self.fingering.append(state) 
165   
166 -    def get_fingerings(self): 
167          " Returns the fingering dictionary of the instrument " 
168          return self.fingerings 
169   
170 -    def set_fingering(self,fingering): 
171          """ Set the fingering of the instrument.  
172   
173              The fingering name must be a key of the fingering dictionary. 
174   
175          """ 
176          if self.fingerings is not None: 
177              if self.fingerings.has_key(fingering): 
178                  self.fingering = self.fingerings[fingering] 
179              else: 
180                  print "Unknown fingering: %s"%fingering 
181          else: 
182              print "Can't set fingering, there is no fingering dictionary" 
183          return self 
184   
185 -    def set_register_key(self,state=1): 
186          " Allow setting the state of the register key " 
187          # this function assumes that the first tonehole of the list is the 
188          # register hole 
189          # this is not always true, and some instruments has many register holes 
190          # a solution must be found 
191          try: 
192              self.fingering[0] = state 
193          except AttributeError: 
194              print "no fingering list available" 
195          return self 
196   
197 -    def get_number_of_sections(self): 
198          " Returns the number of sections the instrument is made off " 
199          return numpy.size(self.sections) 
200   
201 -    def calculate_total_length(self): 
202          " Return the physical length of the instrument. " 
203          L = 0 
204          for section in self.sections: 
205              L = L + section.getLength() 
206          return L 
207   
208 -    def find_toneholes_in_range(self,xi,xf): 
209          """ Return a list of the toneholes located in a specific section of the  
210              instrument  
211          """ 
212          found_toneholes = [] 
213          if self.toneholes is not None: 
214               found_toneholes = [(i,tonehole) for (i,tonehole) in \ 
215                                  enumerate(self.toneholes) if \ 
216                                  xi < tonehole.get_position() <= xf] 
217          return found_toneholes 
218   
219 -def CalculateImpedance(instrument,freqs,T=25.): 
220      """ Calculate the input impedance of the instrument. 
221   
222          This function performs the transfer matrix multiplication. 
223           
224          The frequency vector f must be sufficiently fine for the impedance to 
225          be well represented, otherwise, the CalculateMaxima functions won't be 
226          able to find the maxima. A good vector is: 
227   
228          >>> freqs = numpy.arange(50,5000,5,dtype=float) 
229   
230          @param instrument: an Instrument instance 
231          @param freqs: frequency vector 
232          @param T: the temperature 
233   
234      """ 
235       
236      # make sur f is an array of dimension 1 
237      f = numpy.asarray(freqs).flatten() 
238       
239      # to calculate the input impedance, we must multiply all the  
240      # transmission matrices 
241      # the tonehole transmission matrices must be inserted into the bore sections 
242      # bore sections must therefore be splitted 
243       
244      xi = 0. 
245      TM = create_identity_matrix_array(numpy.size(f)) 
246   
247      #print "Calculating input impedance of the instrument: %s..."%instrument.name 
248      for i,s in enumerate(instrument.sections): 
249          #print 'processing section %d'%i 
250          xf = xi + s.get_length() 
251          found_toneholes = instrument.find_toneholes_in_range(xi, xf) 
252          #print "Number of toneholes in section %d: %d"%(i,len(SubToneHolesList)) 
253   
254          x0 = xi 
255          for (j,tonehole) in found_toneholes: 
256              #print 'Processing tonehole %d'%j 
257              x1 = tonehole.get_position() 
258              TM = multiply_matrix_arrays(TM, s.subTM(f,T,x0-xi,x1-xi)) 
259              try: 
260                  TMh = tonehole.TM(f,s.radius(x1-xi),T,instrument.fingering[j]) 
261              except: 
262                  print "Error processing tone hole no. %d"%j 
263                  raise 
264   
265              TM = multiply_matrix_arrays(TM, TMh) 
266              x0 = x1 # last tonehole become first coordinate of new sub section 
267   
268          # last subsection  
269          # (if there were no subsections, this will be the complete section 
270          TM = multiply_matrix_arrays(TM, s.subTM(f,T,x0-xi,xf-xi)) 
271          xi = xf 
272   
273      TM = multiply_matrix_arrays(TM, instrument.radiation.TM(f,T)) 
274      Z = TM[0]/TM[1] 
275      #print "...done" 
276      return Z[0] 
277   
278 -def CalculateMaxima(Z,inst,freqs,T,number_of_maxima=None): 
279      """ Calculate the exact location of the impedance maxima of an Instrument. 
280   
281          It uses an impedance vector Z, previously calculated with the  
282          CalculateImpedance function, to determine the approximate location 
283          of the maxima and then, uses an optimization algorithm to find their  
284          exact location. 
285   
286          @param Z: impedance vector 
287          @param inst: Instrument instance 
288          @param freqs: frequency vector 
289          @param T: temperature 
290          @param number_of_maxima: optional argument to limit the number of  
291              maxima that are returned 
292          @return: a list of tuples containing frequency and magnitude of each 
293             impedance maxima. 
294      """ 
295   
296      derivative = numpy.diff(numpy.abs(Z)) 
297      signarray = numpy.sign(derivative) 
298      crossing  = numpy.diff(signarray) 
299      indices   = crossing.nonzero() 
300      maximalist = [] 
301       
302      for i in indices[0]: 
303          # do not take into account the minima 
304          if crossing[i] > 0: 
305              continue 
306   
307          if len(maximalist) == number_of_maxima: 
308              break 
309   
310          f0 = freqs[i-1] 
311          f1 = freqs[i+1] 
312   
313          infodict = scipy.optimize.brent( 
314                  lambda f: 1./numpy.abs(CalculateImpedance(inst,f,T)), 
315                  brack=(f0,f1),full_output=True) 
316   
317          fres = infodict[0][0] 
318          zres = 1./infodict[1][0] 
319   
320          # check that it is a maxima 
321          if zres > 1.: 
322              maximalist.append((fres, zres)) 
323   
324      return maximalist 
325