Microwave filters GUI 2.0.3

- Main Page
- Related Pages
- Packages
- Classes
- Files

- File List

libcommonfunc.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3
4 ## @file
5
6 ##
7 # @package libcommonfunc Library of common functions for microwave filters synthesis
8 # @author J.M. Rius, J. Mateu, J.M. Tamayo, C. Collado, A. Padilla and J. O'Callaghan
9 # <p>Dpt. Signal Tehory and Communications, Universitat Politècnica de Catalunya (UPC).
10 # @version 2.0.3
11 # @date June 11, 2013
12 # <p>Acknowledgement: The software was developed in the frame of contract 21398/08/NL/GLC with the European Space Agency (ESA). Technical Offer was Christoph Ernst.
13 # Further features were developed under contract UPC-C7767 with Thales Alenia Space España (TAS-E).
14 # Contributions to the definition of the software functionality and testing have been made by Christoph Ernst, Mónica Martínez Mendoza and other
15 # ESA-ESTEC personnel, and Santiago Sobrino and Luis Roglá from TAS-E.
16 # <p>Copyright: &copy; Universitat Politècnica de Catalunya (UPC) 2009. <br>Contact: lossyfilters@tsc.upc.edu
17 # <p>This program is free software; you can redistribute it and/or modify it under the terms of
18 # the GNU General Public License as published by the Free Software Foundation; either
19 # version 3 of the License, or (at your option) any later version.
20 # <p>This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
21 # without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
22 # See the GNU General Public License for more details.
23 # <p>You should have received a copy of the GNU General Public License along with this program
24 # in file "LICENSE.GPL3"; if not, download it from http://www.gnu.org/licenses/gpl-3.0.html .
25 # <p>Additional permission under GNU GPL version 3 section 7:
26 # <p>If you modify this Program, or any covered work, by linking or combining it with the "Extra Filters Library" (libextrafilters),
27 # the "Lossy Filters Library" (liblossyfilters) or the "License Check Library" (libchecklicense), or modified versions of that libraries,
28 # containing parts covered by the terms described in files LICENSE.LIBEXTRAFILTERS and LICENSE.LIBLOSSYFILTERS,
29 # the licensors of this Program grant you additional permission to convey the resulting work.
30 #
31
32 import re
33 import sys
34 import platform
35 import codecs
36 import os.path
37 import types
38 import csv
39 import numpy as np
40 import scipy as sc
41 import scipy.signal
42 import scipy.optimize as scop
43 from math import log10, pi, atan2, sqrt, floor
44 from PyQt4.QtGui import QMessageBox
45
46 from template import fileTemplate
47
48 # Disable 'Warning: divide by zero encountered in divide' in numpy. Result will be np. inf.
49 # These warning appears tipically when computing Q of resonators for a lossless filter.
50 np.seterr(divide='ignore')
51
52 #Multiprecision tests with mpmath
53 #import mpmath as mp
54 #mp.mp.dps = 100
55
56 #!!!!!!!!!!!!!!!!!!!!!!!!#
57 # Global variables: #
58
59 ##
60 # Dictionary of parameter keywords and variable names.
61 # This dictionary allows the developer to change a variable name in the code without having to change the corresponding keyword in all the existing parameter files.
62 parameterDict = { 'outName': 'outName', 'outDir': 'outDir', 'outDirName': 'outDirName', 'N': 'N', 'f0': 'f0', 'BW': 'BW',
63 'filterTransferFunc': 'filterTransferFunc', 'chebyLr': 'chebyLr', 'qeZero': 'qeZero' , 'qeLr': 'qeLr',
64 'transImagZeros': 'transImagZeros', 'transComplexZeros': 'transComplexZeros', 'genChebyLr': 'genChebyLr',
65 'genChebyLP': 'genChebyLP', 'LPalgor' : 'LPalgor', 'LPmaxIter' : 'LPmaxIter', 'LPNsamples' : 'LPNsamples',
66 'LPcomplexZeros': 'LPcomplexZeros', 'LPrealZeros': 'LPrealZeros', 'LPfracBW': 'LPfracBW', 'LPripple': 'LPripple',
67 'synthTech': 'synthTech', 'zerosArrang': 'zerosArrang', 'adapPredis': 'adapPredis', 'Qef': 'Qef', 'wFunc': 'wFunc',
68 'finiteQo': 'finiteQo', 'Qo': 'Qo', 'QoPredis': 'QoPredis',
69 'nuqK': 'nuqK', 'nuqK11c1': 'nuqK11c1', 'nuqK22c1': 'nuqK22c1', 'nuqK21c1': 'nuqK21c1', 'nuqK21c2': 'nuqK21c2', 'nuqK21c3': 'nuqK21c3',
70 'nuqTech': 'nuqTech', 'nuqDelta': 'nuqDelta',
71 'maxFreq': 'maxFreq', 'minFreq': 'minFreq', 'numFreq': 'numFreq' }
72
73 ##
74 # Number of significant digits to save in files
75 saveSignificantDigits = 14
76
77 ##
78 # Number of significant digits to save in energy and power results.
79 # Usually less than saveSignificantDigits, because not much precision is needed.
80 saveSignificantDigitsEnergy = 4
81
82 ##
83 # Module global variable equal either to None or to the GUI QMainWindow class instance, containing the application main window.
84 # It is equal to None (default) when the calling program is the console interface (mwfiltersgui.py --nogui).
85 # When creating the main window, the GUI executes the libcommonfunc.setGlobalVariablesLibCommonFunc() function in order to set this variable.
86 # The GUI methods do not use this variable. They access the mainWindow class through ParamEditDlg.mainWindow and MainWindow 'self' attributes.
87 mainWindow = None
88
89 ##
90 # Module global variable equal to True if the GUI has been called with the -v flag command line option. Otherwise it is False.
91 verbose = False
92
93 ##
94 # Module global variable equal to True if the GUI has been called with the --nogui flag command line option. Otherwise it is False.
95 nogui = False
96
97 ##
98 # String containing console output generated before the GUI main window exists
99 preGuiOutput = ''
100
101 # End global variables #
102 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
103
104
105 ##
106 #
107 # Parse error catch.
108 #
109 class parseError(Exception):
110 pass
111
112
113 ##
114 #
115 # Synthesis error catch.
116 #
117 class synthError(Exception):
118 pass
119
120
121 ##
122 #
123 # License error catch.
124 #
125 class licenseError(Exception):
126 pass
127
128
129 #!!!!!!!!!!!!!!!!!!!!!!!!#
130 # Global functions: #
131
132 _cfunc_complex = np.vectorize(complex)
133
134 #def mp2np(M):
135 # """
136 # Function to convert an mpmath complex matrix into a numpy array.
137 # @param M = mpmath complex matrix
138 # @returns npM = numpy array
139 # """
140 # npM = np.array( M.tolist() )
141 # return _cfunc_complex(npM.ravel()).reshape(npM.shape)
142 #
143 #
144 #def np2mp(M):
145 # """
146 # Function to convert a numpy complex array or matrix into a mpmath matrix.
147 # @param M = numpy complex array
148 # @returns npM = mpmath matrix
149 # """
150 # return mp.matrix( M.tolist() )
151
152
153 ##
154 #
155 # Function to allow the GUI set values to the global variables of libcomonfunc module.
156 # This function is used in order to avoid circular imports.
157 # @param variablesDict: Variable keyword argument list. Python creates a dictionary.
158 #
159 def setGlobalVariablesLibCommonFunc(**variablesDict):
160 for key, value in variablesDict.items():
161 globals()[key] = value
162
163
164 ##
165 #
166 # Returns the value of the global verbose variable.
167 # Useful to get the value from other modules that were loaded before the value of verbose was updated
168 #
169 def readVerbose():
170 return verbose
171
172
173 ##
174 #
175 # Retrieves software package name and version from string.
176 # @param st = String to parse for package name and version. It must be the docstring (__doc__) of the calling module.
177 # @return pkgname = Package name.
178 # @return version = Package version.
179 #
180 def pkgNameVersion(st):
181 pkgname = re.search('@package (.*)', st).group(1)
182 version = re.search('@version (.*)', st).group(1)
183 return(pkgname, version)
184
185
186 ##
187 #
188 # Prints on standard output if console application.
189 # or on log window if GUI application.
190 # @param st = String.
191 #
192 def myPrint(st):
193
194 # Remove html tags from string:
195 st = re.sub('<[^>]*?>','', st)
196
197 if mainWindow is not None:
198 # In windows platform, remove '\r' end-of-line because Qt expects only \n
199 if platform.system() in [ 'Windows', 'Microsoft', 'Vista' ]: st = st.replace('\r\n', '\n')
200
201 # Another way to do the same:
202 # if sys.platform.find('win') > -1: st = filter(lambda c: c!='\r', st )
203
204 mainWindow.log_listWidget.addItem(unicode(st)) # addItem() method expects a string, not a Parameter class instance
205 mainWindow.log_listWidget.setCurrentRow(mainWindow.log_listWidget.count()-1) # Set log_listWidget current row to last one
206
207 else:
208 # Check if console encoding is able to print non-ascii characters
209 if hasattr(sys.stdout, 'encoding') and sys.stdout.encoding in [None, 'ascii']: # If it is not
210 # st = filter(lambda c: c<128, st ) # Remove non-ascii characters
211 st = ''.join(map(lambda c: c if ord(c)<128 else '?', st )) # Replace non-ascii characters by '?'
212 print st
213
214
215 ##
216 #
217 # Prints application header info.
218 # The information shown is read from the file docstring.
219 # @param applicationName = Name of application, usually os.path.basename(sys.argv[0]).
220 # @param st = String to parse for package name and version. It must be the docstring (__doc__) of the calling module.
221 # @param importModule = Flag, set to True when this function is called to acknowledge successful module import.
222 #
223 def printHeader(applicationName, st, importModule = False):
224
225 if importModule:
226 st = 'Successfully imported: %s; version %s\n' % pkgNameVersion(st)
227 else:
228 st = u"\n%s\n%s" % (applicationName, st)
229 st = st.replace('@package', 'Package')
230 st = st.replace('@author', 'Authors:')
231 st = st.replace('@date', 'Date:')
232 st = st.replace('@version', 'Version:')
233 st = st.replace('<p>', '\n')
234 st = st.replace('<br>', '\n')
235 st = st.replace('&copy;', '(c)')
236 st = """%s \nPython %s on %s \n""" % (st, platform.python_version(), platform.system())
237
238 myPrint(st)
239
240 # Append st to preGuiOutput if the GUI will be opened later
241 if mainWindow is None and not nogui:
242 global preGuiOutput
243 preGuiOutput = preGuiOutput + st
244
245
246 ##
247 #
248 # Read parameters file and return a Parameter class instance.
249 # @param fileName = Parameter file name.
250 # @return Parameter class instance.
251 #
252 def readParamFile(fileName):
253 P = None
254
255
256 try:
257 if fileName:
258 for codec in ['utf-8', 'iso8859_15']:
259 try:
260 myPrint('Openning file: %s' % os.path.basename(fileName))
261 file = codecs.open(fileName, 'r', codec)
262 except IOError, err:
263 raise parseError, 'I/O error: %s' % err
264
265 try:
266 stList = file.readlines()
267 break
268 except UnicodeDecodeError, err:
269 #raise parseError, "Cannot open file with encoding '%s'<p>%s" % (codec, err)
270 myPrint("Cannot open file with encoding '%s'" % (codec))
271
272 myPrint("File sucessfully opened with encoding '%s'" % (codec))
273
274 file.close()
275
276 myPrint('Reading parameters from file: %s' % os.path.basename(fileName))
277 P = Parameters(stList)
278
279 except parseError, err:
280 criticalErrorMsg( 'ERROR in parameter file.<p>' + unicode(err) )
281
282 return P
283
284
285 ##
286 #
287 # Show critical error message. Prints error and exists in console mode. Displays QMessageBox in GUI mode.
288 # @param st = Error string. It can contain <p> and <br> html tags, that are discarded in console mode.
289 #
290 def criticalErrorMsg(st):
291 if mainWindow is not None:
292 QMessageBox.critical(mainWindow, "%s" % mainWindow.applicationName, st)
293
294 ##@todo: See if MainWindow.cleanup() method can be used here.
295 # Close the [S] parameters plot and set all results to None
296 if mainWindow.SPlotMagPhase is not None and not mainWindow.SPlotMagPhase.closed:
297 mainWindow.SPlotMagPhase.close()
298 mainWindow.SPlotMagPhase = None
299 if mainWindow.SPlotError is not None and not mainWindow.SPlotError.closed:
300 mainWindow.SPlotError.close()
301 mainWindow.SPlotError = None
302 CP = CM = SP = None
303 else:
304 st = st.replace('<p>', '\n')
305 st = st.replace('<br>', ' ')
306 st = re.sub('<[^>]*?>','', st)
307 print st
308 sys.exit(2)
309
310
311 ##
312 #
313 # Show warning message. Prints warning in console mode. Displays QMessageBox in GUI mode.
314 # @param st = Warning string. It can contain <p> and <br> html tags, that are discarded in console mode.
315 #
316 def warningMsg(st):
317 if mainWindow is not None:
318 QMessageBox.warning(mainWindow, "%s" % mainWindow.applicationName, st)
319 else:
320 st = st.replace('<p>', '\n')
321 st = st.replace('<br>', ' ')
322 st = st.replace('<b>', '* ')
323 st = st.replace('</b>', ' *')
324 print st
325
326
327 ##
328 #
329 # Show information message. Prints message in console mode. Displays QMessageBox in GUI mode.
330 # @param st = Information string. It can contain <p> and <br> html tags, that are discarded in console mode.
331 #
332 def informationMsg(st):
333 if mainWindow is not None:
334 QMessageBox.information(mainWindow, "%s" % mainWindow.applicationName, st)
335 else:
336 st = st.replace('<p>', '\n')
337 st = st.replace('<br>', ' ')
338 print st
339
340
341 ##
342 #
343 # Show dialog to ask question.
344 # @param st = Information string. It can contain <p> and <br> html tags, that are discarded in console mode.
345 # @param buttons = QMessageBox buttons to show. (e.g.: QMessageBox.Yes | QMessageBox.No)
346 # @return reply = button clicked ((e.g.: QMessageBox.Yes)
347 #
348 def askQuestion(st, buttons):
349 if mainWindow is not None:
350 return QMessageBox.question(mainWindow, "%s" % mainWindow.applicationName, st, buttons)
351 else:
352 pass
353
354
355 ##
356 #
357 # Nicely convert list of floats to string, rounding to 'prec' significant digits
358 # @param listX = List of floats.
359 # @param prec = Number of significant digits (int). Default 3.
360 # @return st = String representation of x.
361 #
362 def listStr(listX, prec=3):
363 st = '['
364 for x in listX: st += '%.*g, ' % (prec, x)
365 return st[0:len(st)-2] + ']'
366
367
368 ##
369 #
370 # Nicely convert complex number to string, rounding real and imaginary parts to 'prec' significant digits
371 # @param x = Complex number.
372 # @param prec = Number of significant digits (int).
373 # @param format = None | 'DB' | 'MA' | 'RI' (string). Default is None.
374 # @param minWidth = Minimum field field (int). Default is 0.
375 # @param eng = Use engineering notation in formats None and 'RI'. Default is False.
376 # Formatted strings of the form 0.x are allowed for exponents 0 and 9.
377 # @param listFormat = Flag to return a list of strings instead of a string. The list contains [ Magnitude, Angle] or [Real, imag]. Useful to save in .csv format.
378 #
379 # The format is borrowed from TouchStone s2p [S] parameter file definition.
380 # None: Complex number, like 2.4-j5.3.
381 # DB: Decibel Magnitude/Degree Angle.
382 # MA: Linear Magnitude/Degree Angle (polar form).
383 # RI: Real Part/Imaginary Part (rectangular form).
384 # @return st = String representation of x.
385 #
386 def complexStr(x, prec, format=None, minWidth=0, eng = False, listFormat = False):
387
388 re, im = np.real(x), np.imag(x)
389
390 # Engineering notation
391 if eng:
392 if re != np.inf:
393 reExp = int(3*floor(log10(abs(re))/3)) if re != 0 else 0
394 re = re / pow(10, reExp)
395 if reExp in [-3, 6]:
396 if re >= 100:
397 re /= 1e3
398 reExp += 3
399 stReExp = 'e' + str(reExp) if reExp != 0 else ''
400 else: stReExp = ''
401
402 if im != np.inf:
403 imExp = int(3*floor(log10(abs(im))/3)) if im != 0 else 0
404 im = im / pow(10, imExp)
405 if imExp in [-3, 6]:
406 if im >= 100:
407 im /= 1e3
408 imExp += 3
409 stImExp = 'e' + str(imExp) if imExp != 0 else ''
410 else: stImExp = ''
411 else:
412 stReExp, stImExp = '', ''
413
414 reim = abs(re/im) if abs(im) > 0 else np.inf
415 imre = abs(im/re) if abs(re) > 0 else np.inf
416
417 if format == None:
418 st1 = '%+.*g%s' % (prec, re, stReExp) if min(abs(re), reim) > 10**-prec else ''
419 st2 = '%+.*g%sj' % (prec, im, stImExp) if min(abs(im), imre) > 10**-prec else ''
420 elif format == 'RI':
421 st1 = '%-+*.*g%s' % (minWidth, prec, re, stReExp) if min(abs(re), reim) > 10**-prec else '%-*s' % (minWidth, '0')
422 st2 = ' %-+*.*g%s' % (minWidth, prec, im, stImExp) if min(abs(im), imre) > 10**-prec else ' %-*s' % (minWidth, '0')
423 st2 = st2.replace('+', ' ')
424 elif format == 'MA':
425 mag, angle = abs(x), (180/pi)*atan2(np.imag(x), np.real(x))
426 st1 = '%-+*.*g' % (minWidth, prec, mag) if abs(mag) > 10**-prec else '%-*s' % (minWidth, '0')
427 st2 = ' %-+*.*g' % (minWidth, prec, angle) if abs(angle) > 10**-prec else ' %-*s' % (minWidth, '0')
428 elif format == 'DB':
429 if x==0: # Set zero to -1000 dB
430 mag = -1000
431 angle = 0
432 else:
433 mag, angle = 20*log10(abs(x)), (180/pi)*atan2(np.imag(x), np.real(x))
434 st1 = '%-+*.*g' % (minWidth, prec, mag) if abs(mag) > 10**-prec else '%-*s' % (minWidth, '0')
435 st2 = ' %-+*.*g' % (minWidth, prec, angle) if abs(angle) > 10**-prec else ' %-*s' % (minWidth, '0')
436
437 if listFormat == False:
438 st = '%-*s' % (minWidth, st1 + st2)
439 if len(st)==0: st = '0'
440 if st[0]=='+': st = ' ' + st[1:len(st)]
441 else:
442 st = [st1, st2]
443
444 return st
445
446 # End global functions
447 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
448
449
450 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
451 # Begin class: Parameters
452
453 ##
454 #
455 # Filter, synthesis and sweep parameters.
456 #
457 class Parameters(object):
458
459 ##
460 #
461 # Parse int, float, complex or string parameters from string stValue.
462 # @param key = string containing parameter keyword.
463 # @param type = string containing parameter type: 'int', 'float', 'complex' or 'string'. 'Float' type allows a value of Inf, inf or INF.
464 # @param stValue = string containing parameter or a list of comma-separated parameters.
465 # @n Complex number format is 1+2j.
466 # @return parameter value.
467 #
468 def readParam(self, key, type, stValue):
469
470 # Remove trailing newline
471 if stValue.endswith('\n'): stValue = stValue[:-1]
472
473 # Check if parametr is empty
474 if stValue.lstrip() == '': return None
475
476 # Split the parameter values list
477 stValue = stValue.split(',')
478
479 # Parse stValue
480 try:
481 if type == 'int': n = [ int(elem) for elem in stValue ]
482 elif type == 'float': n = [ np.inf if elem.lower().lstrip().rstrip()=='inf' else float(elem) for elem in stValue ]
483 elif type == 'complex': n = [ complex(elem) for elem in stValue ]
484 elif type == 'string': n = [ unicode(elem.lstrip().rstrip(), 'utf-8') for elem in stValue ]
485 else: raise parseError, "Bad [%s] parameter type: '%s' " % (key, type)
486 except ValueError, err:
487 raise parseError, "Bad [%s] parameter format: '%s' " % (key, err)
488
489 # Convert one-element list to element type
490 if len(n) == 1: n = n[0]
491
492 return n
493
494
495 ##
496 #
497 # Initialize parameters class with list of strings.
498 # Each string is a line of the parameter file.
499 # For backwards compatibility of parameter files, missing fields are set to None.
500 # @param stList = List of strings. Each element of list is a line of input file.
501 # @n If stList is equal to None, an empty instance of Parameters class is created.
502 # (useful to create a new parameter file and for the copy funcion).
503 # @return Parameter class instance.
504 #
505 def __init__(self, stList):
506 if stList is not None: # Create Parameter class instance from parameter information in stList
507
508 # Extract (keyword, type, value) from lines containing parameter values
509 validLines = [ st for st in stList if st.find('=') >=0 and not st.startswith('#') ]
510
511 pat = re.compile('^.*\[(.+)\].*\{(.+)\}.*=(.+)$')
512 matchList = [ pat.match(line) for line in validLines ]
513
514 if matchList[0] is None: # Chebash file format
515 myPrint("Reading Chebash file format")
516 pat = re.compile('^(.+)=(.+)$')
517 matchList = [ pat.match(line) for line in validLines ]
518 chebashFormat = True
519
520 else: # LossyFilters file format
521 myPrint("Reading LossyFilters file format")
522 chebashFormat = False
523
524 # List of tuples containing (keyword, type, value) for each parameter
525 paramList = [ match.groups(1) for match in matchList if match is not None]
526
527 if chebashFormat: # Read paramList in Chebash fornat
528 self.readChebash(paramList)
529 else: # Read paramList in Lossyfilters fornat
530 try:
531 for (key, type, value) in paramList:
532 # Double backslashes have been removed by parameter matching, and must be restored. Otherwise, a utf8 decoding error will result from \nnn secuences
533 # It is necessary also to strip \r when reading windows files.
534 value = value.replace('\\', '\\\\').strip()
535 command = "self.%s = self.readParam('%s', '%s', '%s')" % ( parameterDict[key], key, type, value )
536 exec command
537 except KeyError: raise parseError, u"Invalid keyword: '%s'" % key
538
539 # For backwards compatibility of parameter files, set missing fields to None.
540 for field in parameterDict.values():
541 if field not in dir(self): exec "self.%s = None" % field
542
543 #For backwards compatibility with Parameter files prior to implementation of Linear Phase optimization parameters
544 if self.LPalgor is None: self.LPalgor = 'SLSQP'
545 if self.LPmaxIter is None: self.LPmaxIter = 1000
546 if self.LPNsamples is None: self.LPNsamples = 100
547
548 # For backwards compatibility with Parameter files prior to implementation of asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11):
549 if self.nuqK11c1 is None: self.nuqK11c1 = self.nuqK
550 if self.nuqK22c1 is None: self.nuqK22c1 = self.nuqK
551 if self.nuqK21c1 is None: self.nuqK21c1 = self.nuqK
552 if self.nuqK21c2 is None: self.nuqK21c2 = self.nuqK
553 if self.nuqK21c3 is None: self.nuqK21c3 = self.nuqK
554
555 # Check if parameters have invalid values and make some basic parameter processing
556 self.validate()
557
558 ##
559 # Flag that indicates if the class instance contains information.
560 # empty==False or empty==True
561 self.empty = False
562
563 else: # Create empty Parameter class instance with all fields set to None
564 for key in parameterDict.values():
565 command = "self.%s = None" % key
566 exec command
567 self.empty = True
568 # Set lists to empty list:
569 self.transImagZeros = []
570 self.transComplexZeros = []
571 self.LPcomplexZeros = []
572 self.LPrealZeros = []
573 self.wFunc = []
574
575 # Set values no to crash the GUI
576 self.N = 2
577 self.zerosArrang = 3
578
579 ##
580 # Path of output files, including directory and file name
581 self.outDirName = None
582
583
584 ##
585 #
586 # Translate a list of Chebash parameters into a list of equivalent lossyfilters parameters
587 # @param paramList = List of tuples (Chebash_Parameter_Name, value) where value is a string
588 # @return paramList = List of tuples (LossyFilters_Parameter_Name, type, value) where type and value are strings
589 #
590 def readChebash(self, paramList):
591
592 ## def normRealFreq(f):
593 ## """
594 ## Auxiliary function
595 ## """
596 ## return (f/self.f0 - self.f0/f)*self.f0/self.BW
597
598
599 ## for (keyChebash, stValue) in paramList:
600
601 # We have to assume that Chebash parameters appear in the same line order.
602 # We cannot use keywords to find the parameters, since there are repeated keywords (e.g. 'real', 'imag')
603 if paramList[0][0] == 'n': self.N = int(paramList[0][1])
604 else: raise parseError, u"parameter 'n' missing or misplaced in Chebash file"
605
606 if paramList[1][0] == 'fo': self.f0 = float(paramList[1][1]) * 1e6
607 else: raise parseError, u"parameter 'fo' missing or misplaced in Chebash file"
608
609 if paramList[2][0] == 'Qo': self.Qo = float(paramList[2][1])
610 else: raise parseError, u"parameter 'Qo' missing or misplaced in Chebash file"
611
612 if paramList[3][0] == 'wp': self.BW = float(paramList[3][1]) * 1e6
613 else: raise parseError, u"parameter 'wp' missing or misplaced in Chebash file"
614
615 if paramList[4][0] == 'ro': self.genChebyLr = float(paramList[4][1])
616 else: raise parseError, u"parameter 'ro' missing or misplaced in Chebash file"
617
618 if paramList[5][0] == 'puntos': self.numFreq = int(paramList[5][1])
619 else: raise parseError, u"parameter 'puntos' missing or misplaced in Chebash file"
620
621 if paramList[6][0] == 'n_ceros': n_ceros = int(paramList[6][1])
622 else: raise parseError, u"parameter 'n_ceros' missing or misplaced in Chebash file"
623
624 if n_ceros > 0:
625 if paramList[7][0] == 'ceros': self.transImagZeros = [ self.f0 + float(st)*1e6 for st in paramList[7][1].split(',') ]
626 else: raise parseError, u"parameter 'ceros' missing or misplaced in Chebash file"
627 if len(self.transImagZeros) != n_ceros: raise parseError, u"Incorrect number of transmission zeros in Chebash file"
628 else:
629 self.transImagZeros = []
630
631 if paramList[8][0] == 'n_ceros_eq': n_ceros_eq = int(paramList[8][1])
632 else: raise parseError, u"parameter 'n_ceros_eq' missing or misplaced in Chebash file"
633
634 if n_ceros_eq > 0:
635 ## if paramList[9][0] == 'real': realEqZs = [ normRealFreq(self.f0 + float(st)*1e6) for st in paramList[9][1].split(',') ]
636 if paramList[9][0] == 'real': realEqZs = [ float(st) for st in paramList[9][1].split(',') ]
637 else: raise parseError, u"parameter 'real' (for equalization ceros) missing or misplaced in Chebash file"
638
639 if paramList[10][0] == 'imag': imagEqZs = [ self.f0 + float(st)*1e6 for st in paramList[10][1].split(',') ]
640 else: raise parseError, u"parameter 'imag' (for equalization ceros) missing or misplaced in Chebash file"
641
642 self.transComplexZeros = [ complex(val[0],val[1]) for val in zip(realEqZs, imagEqZs) ]
643 if len(self.transComplexZeros) != n_ceros_eq: raise parseError, u"Incorrect number of equalization zeros in Chebash file"
644 else:
645 self.transComplexZeros = []
646
647 # Remove complex transmission zeros with symmetric real part
648 if n_ceros_eq > 0:
649 nonSymmZeros = [ ] # List of zeros removing symmetric in the real part
650 for TZ in self.transComplexZeros:
651 if -TZ.conjugate() not in nonSymmZeros: nonSymmZeros.append(TZ)
652 self.transComplexZeros = nonSymmZeros
653
654 # Default values for variables not in Chebash files
655 self.outName = 'Chebash'
656 self.outDir = '.'
657 self.filterTransferFunc = 'Generalized Chebyshev'
658 self.genChebyLP = 0
659 self.synthTech = 'Minimum Insertion Loss'
660 self.finiteQo = 1
661 self.nuqTech = 'kS21+kS11'
662 self.minFreq = self.f0 - 3*self.BW
663 self.maxFreq = self.f0 + 3*self.BW
664
665 ## ##
666 ## # List of Chebash parameters to ignore
667 ## chebashIgnoreList = ['n_ceros', 'n_frecpropias_eq', 'real', 'imag', 'tipo_transformacion', 'graf[1]', 'graf[2]', 'graf[3]', 'graf[5]', 'graf[6]', 'graf[7]',
668 ## 'graf[8]', 'graf[9]', 'graf[10]', 'graf[11]', 'graf[12]', 'graf[13]', 'graf[14]', 'graf[15]', 'graf[16]', 'graf[17]', 'graf[18]']
669 ##
670 ## ##
671 ## # List of Chebash Mask parameters not to read
672 ## chebashMaskList = [ 'n_at', 'f_at', 'A_at', 't_at', 'offsetmba', 't_offsetmba', 'n_perd', 'f_abp', 'a_abp', 't_abp', 'offsetmbp', 't_offsetmbp',
673 ## 'nr_retardo', 'f_retardo', 'r_retardo', 't_retardo', 'offsetrg', 't_offsetrg', 'nptos_pendht', 'f_pendht', 'd_pendht', 't_pendht', 'offsetpendht', 't_offsetpendht',
674 ## 'nptos_pendre', 'f_pendre', 'd_pendre', 't_pendre', 'offsetpendre', 't_offsetpendre', 'nptos_gdst', 'f_gdst', 'v_gdst', 't_gdst', 'offsetgdst', 't_offsetgdst',
675 ## 'nptos_gdsp', 'f_gdsp', 'v_gdsp', 't_gdsp', 'offsetgdsp', 't_offsetgdsp', 'n_frecpropias_eq', 'variacion_temp', 'margen_temp', 'modificar_mascaras', 'margen_gdsp']
676
677 ## processingComplex = False # Will be True when we start processing a 3-line complex parameter
678 ## for (keyChebash, stValue) in paramList:
679 ##
680 ## if not processingComplex: # int, float, string or 1st line of a complex parameter
681 ## try:
682 ## (keyLossyFilters, type, scale) = chebashTranslateDict[keyChebash]
683 ## except KeyError, err:
684 ## if keyChebash not in chebashIgnoreList + chebashMaskList: myPrint('Unknown parameter key in Chebash file: %s' % err)
685 ## continue
686 ##
687 ## if type in ['int', 'float']:
688 ##
689 ## # First process the special cases
690 ## if keyChebash == 'graf[4]':
691 ## fmin, fmax = tuple( stValue.split(',')[0:2] )
692 ## paramListLF.append( ('minFreq', 'float', str(float(fmin)*scale) ) )
693 ## paramListLF.append( ('maxFreq', 'float', str(float(fmax)*scale) ) )
694 ## continue
695 ## else:
696 ## if stValue != 'X': # In Chebash, 'X' indicates an empty parameter
697 ## exec("value = [ scale* %s(n) for n in stValue.split(',') ]" % type)
698 ## if len(value) == 1: # Not a list
699 ## value = str(value[0])
700 ## else: # List
701 ## value = str(value).replace('[', '').replace(']', '')
702 ## else: # Empty parameter
703 ## value = ' '
704 ##
705 ## elif type == 'string':
706 ## value = stValue
707 ## elif type == 'complex':
708 ## processingComplex = True
709 ## NComplex = int(stValue)
710 ## else:
711 ## raise parseError, 'Invalid data type'
712 ##
713 ## else: # 2nd and 3rd lines of a complex parameter
714 ## if keyChebash == 'real':
715 ## if NComplex > 0:
716 ## complexRe = [ float(n)*scale for n in stValue.split(',') ]
717 ## elif keyChebash == 'imag':
718 ## if NComplex > 0:
719 ## complexIm = [ float(n)*scale for n in stValue.split(',') ]
720 ## value = str([ complex(nReal, nImag) for nReal, nImag in zip(complexRe, complexIm) ]).replace('[', '').replace(']', '').replace('(', '').replace(')', '')
721 ## else:
722 ## value = ' '
723 ## processingComplex = False
724 ##
725 ## # Add lossyfilters parameter, for normal or for special case
726 ## paramListLF.append((keyLossyFilters, type, value))
727 ##
728 ## return paramListLF + chebashDefaultList
729
730
731 ##
732 #
733 # Function to check if parameters have invalid values and make some basic parameter processing.
734 #
735 def validate(self):
736
737 if not isinstance(self.outName, types.StringTypes) or len(self.outName) == 0:
738 raise parseError, u"Name of results output files is not a valid string"
739
740 # DO NOT CHECK OUTDIR KNOW. CKECK IT AFTER COMPUTATION, JUST BEFORE SAVING
741 # THIS WILL ALLOW MODIFICATION OF OUTDIR USING THE GUI
742 #
743 # if self.outDir is not None:
744 # # Normalize and expand output directory name, just in case
745 # self.outDir = os.path.normpath(self.outDir)
746 # self.outDir = os.path.expanduser(self.outDir)
747 # if not isinstance(self.outDir, types.StringTypes) or not os.path.isdir(self.outDir):
748 # raise parseError, u"Output directory parameter '%s' is not an existing directory name" % self.outDir
749 #
750 # # Create output name using outDir and outName
751 # if self.outDir is not None: self.outDirName = os.path.join(self.outDir, self.outName)
752 # else: self.outDirName = self.outName
753 #
754 # # Normalize output path just in case
755 # self.outDirName = os.path.normpath(self.outDirName)
756
757 if self.filterTransferFunc.title() not in [ u'Butterworth', u'Chebyshev', u'Quasieliptic', u'Generalized Chebyshev' ]:
758 raise parseError, u"Invalid filter topology type '%s'" % self.filterTransferFunc
759
760 if self.synthTech.title() not in [ u'Minimum Insertion Loss', u'Predistortion', u'Prescribed Insertion Loss' ]:
761 raise parseError, u"Invalid synthesis technique '%s'" % self.synthTech
762
763 if self.nuqTech not in [ u'kS21+kS11', u'kS21', u'kS21+pole/zero' ]:
764 raise parseError, u"Invalid Prescribed Insertion Loss technique = '%s'" % self.nuqTech
765
766 if self.numFreq <= 0:
767 raise parseError, u"Sweep parameters: Number of samples is not > 0"
768
769 if self.minFreq >= self.maxFreq:
770 raise parseError, u"Sweep parameters: Minimum frequency is not smaller than Maximum frequency"
771
772 if self.minFreq <= 0:
773 raise parseError, u"Minimum frequency must be larger than 0 Hz"
774
775 if self.f0<= 0:
776 raise parseError, u"Central frequency must be larger than 0"
777
778 if self.BW<= 0:
779 raise parseError, u"Bandwidth must be larger than 0"
780
781 if self.numFreq <= 1:
782 raise parseError, u"Number of frequency samples must be larger than 1"
783
784 if self.filterTransferFunc.title() == u'Chebyshev' and self.chebyLr <= 0:
785 raise parseError, u"Chebyshev filter return losses must be larger than 0"
786
787 if self.filterTransferFunc.title() == u'Quasieliptic':
788 if self.qeLr <= 0: raise parseError, u"Quasieliptic filter return losses must be larger than 0"
789 if self.qeZero <= 1: raise parseError, u"Quasieliptic filter transmission zero must be larger than 1"
790
791 # Lists that could be empty in the parameter file
792 if self.transImagZeros is None: self.transImagZeros = [ ]
793 if self.transComplexZeros is None: self.transComplexZeros = [ ]
794 if self.LPcomplexZeros is None: self.LPcomplexZeros = [ ]
795 if self.LPrealZeros is None: self.LPrealZeros = [ ]
796 if self.wFunc is None: self.wFunc = [ ]
797
798 # Be sure that lists of zeros are really lists
799 if type(self.transImagZeros) != types.ListType: self.transImagZeros = [ self.transImagZeros ]
800 if type(self.transComplexZeros) != types.ListType: self.transComplexZeros = [ self.transComplexZeros ]
801 if type(self.LPcomplexZeros) != types.ListType: self.LPcomplexZeros = [ self.LPcomplexZeros ]
802 if type(self.LPrealZeros) != types.ListType: self.LPrealZeros = [ self.LPrealZeros ]
803
804 # Lists that should be a one-element list when it is a number
805 if type(self.LPcomplexZeros).__name__ in ['int', 'float', 'complex']: self.LPcomplexZeros = [ self.LPcomplexZeros ]
806 if type(self.LPrealZeros).__name__ in ['int', 'float', 'complex']: self.LPrealZeros = [ self.LPrealZeros ]
807
808 if self.filterTransferFunc.title() == u'Generalized Chebyshev':
809 if self.genChebyLr <= 0:
810 raise parseError, u"Generalized Chebyshev filter return losses must be larger than 0"
811
812 if self.genChebyLP and any( np.real(self.LPcomplexZeros + self.LPrealZeros)==0 ):
813 raise parseError, 'Linear Phase optimization zeros cannot be purely imaginary'
814
815 if self.genChebyLP and any( np.imag(self.LPcomplexZeros)<=0 ):
816 raise parseError, 'Imaginary part of Linear Phase optimization complex zeros cannot be negative or 0 Hz'
817
818 if len(self.transImagZeros)>0 and any(np.array(self.transImagZeros)<=0):
819 raise parseError, 'Imaginary transmission zeros cannot be negative or 0 Hz'
820
821 if len(self.transComplexZeros)>0 and any(np.imag(self.transComplexZeros) <=0):
822 raise parseError, 'Imaginary part of complex transmission zeros cannot be negative or 0 Hz'
823
824 if self.synthTech.title() == u'Minimum Insertion Loss':
825 if self.finiteQo and self.Qo <=0: raise parseError, u"Q of implementation must be larger than 0"
826
827 if self.synthTech.title() == u'Predistortion':
828 if self.QoPredis <=0: raise parseError, u"Predistortion Q of implementation must be larger than 0"
829 if self.Qef <=0: raise parseError, u"Predistortion Q to emulate must be larger than 0"
830
831 if type(self.zerosArrang).__name__ == 'int' and self.zerosArrang>=1 and self.zerosArrang<=4: pass
832 elif type(self.zerosArrang).__name__ == 'list':
833 if len(self.zerosArrang) != self.N: raise parseError, u"Length of Predistortion factorization method boolean list is different than filter order"
834 if type(self.zerosArrang[0]).__name__ == 'bool': pass
835 elif type(self.zerosArrang[0]).__name__ == 'int':
836 for n in range(self.N): self.zerosArrang[n] = bool(self.zerosArrang[n])
837 else: raise parseError, u"Predistortion factorization method is not a list of Booleans or integers than can be converted to Booleans"
838
839 else: raise parseError, u"Predistortion factorization method is not an integer in the range [1,4] or a list of Booleans"
840
841 if self.adapPredis == 1:
842 if len(self.wFunc) != self.N: raise parseError, u"Length of adaptive predistortion weights list different than filter order"
843
844 if self.synthTech.title() == u'Prescribed Insertion Loss':
845 if self.nuqTech == 'kS21+kS11':
846 if self.nuqK11c1 <=0: raise parseError, u"Prescribed Insertion Loss k21+kS11: k11 must be larger than 0"
847 if self.nuqK22c1 <=0: raise parseError, u"Prescribed Insertion Loss k21+kS11: k22 must be larger than 0"
848 if self.nuqK21c1 <=0: raise parseError, u"Prescribed Insertion Loss k21+kS11: k21 must be larger than 0"
849 if self.nuqTech == 'kS21':
850 if self.nuqK21c2 <=0: raise parseError, u"Prescribed Insertion Loss k21: k21 must be larger than 0"
851 if self.nuqTech == 'kS21+pole/zero':
852 if self.nuqK21c3 <=0: raise parseError, u"Prescribed Insertion Loss kS21+pole/zero: k21 must be larger than 0"
853 if self.nuqDelta <= 0: raise parseError, u"Prescribed Insertion Loss kS21+pole/zero: Delta must be larger than 0"
854 eps = 0.05 # Tolerance in the validation of Delta parameter
855 minDelta = 20 * 10**(-self.nuqK21c3/10) - eps
856 if self.nuqDelta < minDelta: warningMsg("WARNING: Prescribed Insertion Loss kS21+pole/zero:<p>Delta parameter smaller than the<br>minimum recommended value (%.2g)" % minDelta)
857
858
859 ##
860 #
861 # Function to copy instances of the Parameters class (copies the self instance).
862 # @return Pcopy = New instance of Parameters class copied from the self instance.
863 #
864 def copy(self):
865 Pcopy = Parameters(None)
866
867 if self.empty == False: # If the class instance to copy is not empty, copy it. Otherwise just return the empty Pcopy
868 for key in parameterDict.values():
869 command = "Pcopy.%s = self.%s" % (key, key)
870 exec command
871 Pcopy.empty = False
872
873 return Pcopy
874
875
876 ##
877 #
878 # Special method to print parameter class (prints the self instance).
879 # Floating point variables are printed as strings using str() because str() is inmune to roundoff errors in the floating point representation.
880 #
881 # @return st = Human redable printing of the parameters (string).
882 #
883 def __str__(self):
884 if self.empty == True:
885 print ''
886 return
887
888 # Pattern to find parameter keywords
889 pat = re.compile('^.*\[(.+)\].*')
890
891 out = ''
892 for line in fileTemplate.split('\n'):
893 # Substitue parameter values in template
894 if line.find('=') >=0 and not line.startswith('#'):
895 key = pat.match(line).groups(0)[0]
896 value = eval('self.' + parameterDict[key])
897
898 # Round floats and lists of floats (converted to string in order to represent properly the rounded number)
899 if type(value) == types.FloatType or type(value) == types.ComplexType:
900 value = complexStr(value, saveSignificantDigits, eng = True)
901 elif type(value) == types.ListType:
902 value = ', '.join( [ complexStr(val, saveSignificantDigits, eng = True) for val in value ] )
903
904 if value is not None:
905 st = 'Inf' if value==np.inf else unicode(value).strip()
906 else:
907 st = ''
908
909 if type(value) == types.UnicodeType:
910 line = line + filter(lambda c: c not in "[]", st ) # Allow () in strings
911 else:
912 line = line + filter(lambda c: c not in "()[]", st ) # Remove () from complex numbers
913
914 # Write platform-dependant end-of-line
915 if platform.system() in [ 'Windows', 'Microsoft', 'Vista' ]:
916 out = out + line + '\r\n'
917 elif platform.system() in [ 'Linux' ]:
918 out = out + line + '\n'
919 else: # Mac platform
920 out = out + line + '\r'
921
922 return out
923
924 # End class: Parameters #
925 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
926
927
928 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
929 # Begin class: Sparameters
930
931 ##
932 #
933 # [S] parameters class.
934 #
935 class Sparameters(object):
936
937 ##
938 #
939 # Initialize frequency axis attributes and compute S parameters from Characteristic Polynomials.
940 #
941 # @param P = Parameters class instance, containing filter, synthesis and sweep parameters.
942 # @param CP = CharPolys class instance, containing characteristic polynomials and constants.
943 # @param FT = Frequency transform class instance, containing methods to normalize frequency and unnormalize group delay.
944 # @param freq = Frequency axis (Hz). If equal to None, the frequency sweep based on P.minFreq, P.maxFreq and P.numFreq will be used. Default None.
945 # @param checkPassive = If checkPassive is True, a warning will be raised when there are [S] parameters larger than 0 dB. Default True.
946 #
947 # [S] parameters computed from characteristic polynomials.
948 # \f[
949 # S_{11} = \frac{1}{\epsilon_R} \frac{F(s)}{E(s)}
950 # \f]
951 # \f[
952 # S_{21} = S_{12} = \frac{1}{\epsilon} \frac{P(s)}{E(s)}
953 # \f]
954 # \f[
955 # S_{22} = \frac{1}{\epsilon_R} \frac{F_{22}(s)}{E(s)} = (-1)^N \frac{1}{\epsilon_R} \frac{F(s)^*}{E(s)}
956 # \f]
957 # except for asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11), where it is computed as
958 # \f[
959 # S_{22} = k_{22} \left( 1 - \left( \frac{k_{21}}{k_{11}} \right) ^2 \right) + \left( \frac{k_{21}}{k_{11}} \right) ^2 S_{11}
960 # \f]
961 # In this "Prescribed Insertion Loss technique" case 1, the correct lossy values of \f$ S_{11} \f$ and \f$ S_{21} \f$
962 # are obtained because \f$ \epsilon_R \f$ and \f$ \epsilon \f$ have been divided by \f$ k_{11} \f$ and \f$ k_{21} \f$ , respectively,
963 # when processing the polinomials at liblossyfilters.nonUniformQ() function .
964 #
965 def __init__(self, P, CP, FT, freq = None, checkPassive=True):
966 if CP == None: return
967
968 ##
969 # Fractional bandwidth of the FT used to compute this [S] parameters.
970 self.FBW = FT.FBW
971
972 ##
973 # Central frequency of the FT used to compute this [S] parameters.
974 self.f0 = FT.f0
975
976 ##
977 # Frequency transform used to compute this [S] parameters.
978 self.FT = FT
979
980 ##
981 # Frequency axis (Hz)
982 self.freq = None
983
984 if freq is None:
985 self.freq = np.linspace(P.minFreq, P.maxFreq, P.numFreq)
986 else:
987 self.freq = freq
988
989 ##
990 # Normalized \f$ \omega ' \f$ axis
991 self.omega = FT.normFreq(self.freq)
992
993 # Value of quality factor used to compute dissipation factor sigma.
994 # It is infnite in Prescribed Insertion Loss and in Minimum Insertion Loss ideal filter, otherwise it is the value in parameters file.
995 if P.synthTech == u'Minimum Insertion Loss':
996 if P.finiteQo: Qo = P.Qo
997 else: Qo = np.Inf
998 elif P.synthTech == u'Predistortion':
999 Qo = P.QoPredis
1000 elif P.synthTech == u'Prescribed Insertion Loss':
1001 Qo = np.Inf
1002
1003 ##
1004 # Dissipation factor. [Cameron eq. (3.118)] [TN 101.1 eq. (14)].
1005 #
1006 # \f$ \sigma = \frac{f_0}{\Delta f} \frac{1}{Q_0} \f$
1007 #
1008 self.sigma = 1 / (Qo * FT.FBW)
1009
1010 ##
1011 # Complex frequency variable in normalized s-plane.
1012 # \f$ s = \sigma + j \omega ' \f$
1013 self.s = self.sigma + 1j * self.omega
1014
1015 ##
1016 # Parameter \f$ S_{11} \f$ computed from characteristic polynomials.
1017 # \f$ S_{11} = \frac{1}{\epsilon_R} \frac{F(s)}{E(s)} \f$
1018 self.S11 = None
1019
1020 ##
1021 # Parameter \f$ S_{12} \f$ computed from characteristic polynomials.
1022 # \f$ S_{12} = \frac{1}{\epsilon} \frac{P(s)}{E(s)} \f$
1023 self.S12 = None
1024
1025 ##
1026 # Parameter \f$ S_{21} \f$ computed from characteristic polynomials.
1027 # \f$ S_{21} = \frac{1}{\epsilon} \frac{P(s)}{E(s)} \f$
1028 self.S21 = None
1029
1030 ##
1031 # Parameter \f$ S_{22} \f$ computed from characteristic polynomials.
1032 # \f$ S_{22} = \frac{1}{\epsilon_R} \frac{F_{22}(s)}{E(s)} = (-1)^N \frac{1}{\epsilon_R} \frac{F(s)^*}{E(s)} \f$
1033 # [Cameron eq. (6.15) pp. 212]
1034 #
1035 # except for asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11), where it is computed as
1036 # \f$ S_{22} = 1 - \left ( \frac{k_{21}}{k_{11}} \right) ^2 + \frac{k^2_{21}}{k_{11}k_{22}} S_{11_lossless} \f$
1037 #
1038 self.S22 = None
1039
1040 # Compute self.S11,self.S12, self.S21, self.S22
1041 self.fromCharPolys(P, self.s, CP.Ps, CP.Es, CP.Fs, CP.eps, CP.epsR)
1042
1043 if checkPassive:
1044 epstol = 1e-3
1045 S11max, S21max, S22max = max(abs(self.S11)), max(abs(self.S21)), max(abs(self.S22))
1046 if max([S11max, S21max, S22max]) > 1+epstol:
1047 warningMsg("""WARNING: [S] parameters computed from CP are larger than 0 dB at some frequency points<p>
1048 max[S11] = %.2f dB, max[S21] = %.2f dB, max[S22] = %.2f dB
1049 """ % (20*log10(S11max), 20*log10(S21max), 20*log10(S22max)) )
1050 # if any(abs(self.S11) > 1+epstol) or any(abs(self.S21) > 1+epstol) or any(abs(self.S22) > 1+epstol):
1051 # warningMsg("WARNING: [S] parameters computed from CP are larger than 0 dB at some frequency points")
1052
1053 ##
1054 # Group delay for the LowPass prototype
1055 self.groupDelayLP = self.groupDelayfunc(CP.Es, CP.Ps, self.s)
1056 # Computation of group delay as a derivative of the phase of S21
1057 # Fails when sampling S21 near the transmission zeros
1058 #self.groupDelayLP = - np.ediff1d(self.phaseS21) / np.ediff1d(self.omega)
1059 # Normalized \f$ \omega \f$ at the sampling points of the group delay
1060 #self.omegaGD = 0.5*(self.omega[1:] + self.omega[:-1])
1061
1062 ##
1063 # Group delay for the unnormalized filter (sec)
1064 self.groupDelay = FT.unormGroupDelay(self.freq, self.groupDelayLP)
1065 self.groupDelay[ self.groupDelay<0 ] = 0 # CE prefers that the GD goes negative but TAS-E prefers truncation of negative values
1066
1067 ##
1068 # Phase of transfer function (rad)
1069 self.phaseS21 = np.unwrap(np.angle(self.S21))
1070
1071 ##
1072 # Phase of return losses (rad)
1073 self.phaseS11 = np.unwrap(np.angle(self.S11))
1074
1075 ##
1076 # Phase of inverse return losses (rad)
1077 self.phaseS22 = np.unwrap(np.angle(self.S22))
1078
1079 # Best fit to linear phase straight line, fitting to straight line in the whole band
1080 #coef = np.polyfit(self.omega, self.phaseS21, 1)
1081
1082 if freq is None:
1083 # Best fit to linear phase straight line, fitting to straight line in the center of the band in an interval equal to half the bandwidth
1084 NBW = int( P.numFreq * P.BW/(P.maxFreq - P.minFreq) ) # Number of samples corresponding to BW
1085 N1 = int(P.numFreq/2 - NBW/4)-1
1086 N2 = int(P.numFreq/2 + NBW/4)+1
1087 coef = np.polyfit(self.omega[N1:N2], self.phaseS21[N1:N2], 1)
1088
1089 ##
1090 # Deviation from linear phase.
1091 # It is the error between the phase of S21 and the straight line that best fits the phase in the center of the band in an interval equal to half the bandwidth.
1092 self.deviationLP = self.phaseS21 - (self.omega*coef[0] + coef[1])
1093
1094 ##
1095 # Parameter \f$ S_{11} \f$ computed from a coupling matrix.
1096 self.S11M = None
1097
1098 ##
1099 # Parameter \f$ S_{21} \f$ computed from a coupling matrix.
1100 self.S21M = None
1101
1102 ##
1103 # Parameter \f$ S_{22} \f$ computed from a coupling matrix.
1104 self.S22M = None
1105
1106 ##
1107 # Power dissipated at resonators.
1108 # 2-D numpy array of shape (Nf, Nr), where Nf is the number of frequencies and Nr the number of resonators.
1109 self.powerRes = None
1110
1111 ##
1112 # Energy stored at resonators.
1113 # 2-D numpy array of shape (Nf, Nr), where Nf is the number of frequencies and Nr the number of resonators.
1114 self.energyRes = None
1115
1116 ##
1117 # Power dissipated at resistive couplings.
1118 # Tuple ( (m, n), Gmn, powerDirect, powerReverse ) where power is a 1-D numpy array with one element for each frequency.
1119 # Only entries n>m are stored.
1120 self.powerCoup = None
1121
1122 ##
1123 # Energy stored at magnetic or electric couplings.
1124 # Tuple ( (m, n), Bmn, energyDirect, energyReverse ) where energy is a 1-D numpy array with one element for each frequency.
1125 # Only entries n>m are stored.
1126 self.energyCoup = None
1127
1128 ##
1129 # Power dissipated at resonators, with reverse excitation.
1130 # 2-D numpy array of shape (Nf, Nr), where Nf is the number of frequencies and Nr the number of resonators.
1131 self.powerResRev = None
1132
1133 ##
1134 # Energy stored at resonators, with reverse excitation.
1135 # 2-D numpy array of shape (Nf, Nr), where Nf is the number of frequencies and Nr the number of resonators.
1136 self.energyResRev = None
1137
1138
1139 ##
1140 #
1141 # Compute S-Parameters (\f$S_{11}\f$, \f$S_{12}\f$, \f$S_{21}\f$, and \f$S_{22}\f$) from the characteristic polynomials.
1142 # Evaluation points are taken from attribute Sparameters.s .
1143 # This method will be called by the BW and TZs optimization routines.
1144 # Output in atributes Sparameters.S11, Sparameters.S12, Sparameters.S21 and Sparameters.S22 .
1145 #
1146 # @param P = Parameters class instance, containing filter, synthesis and sweep parameters.
1147 # @param evalS = Vector with the normalized complex frequencies (s) where we want to compute the group delay.
1148 # @param Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$.
1149 # @param Es = Characteristic polynomial of the denominator of \f$S_{11}\f$ and \f$S_{21}\f$.
1150 # @param Fs = Characteristic polynomial of the numerator of \f$S_{11}\f$.
1151 # @param eps = Characteristic constant \f$ \epsilon \f$ in the denominator of \f$ S_{12} \f$ and \f$ S_{21} \f$
1152 # @param epsR = Characteristic constant \f$ \epsilon_R \f$ in the denominator of \f$ S_{11} \f$
1153 #
1154 def fromCharPolys(self, P, evalS, Ps, Es, Fs, eps, epsR):
1155
1156 # Evaluate polynomials for s: The result is of type numpy.array
1157 aFs = Fs(evalS)
1158 aPs = Ps(evalS)
1159 aEs = Es(evalS)
1160
1161 self.S11 = (1/epsR) * aFs / aEs
1162 self.S12 = (1/eps) * aPs / aEs
1163 self.S21 = self.S12
1164
1165 if P.synthTech.title() == u'Prescribed Insertion Loss' and P.nuqTech == u'kS21+kS11':
1166 k11 = 10**(-P.nuqK11c1/20.0) # Insertion lossess
1167 k22 = 10**(-P.nuqK22c1/20.0) # Insertion lossess
1168 k21 = 10**(-P.nuqK21c1/20.0) # Insertion lossess
1169 self.S22 = k22 * ( 1 - (k21/k11)**2 ) + ( k21/k11)**2 * (self.S11)
1170 else:
1171 self.S22 = (1/epsR) * (-1)**P.N * aFs.conj() / aEs
1172 # Do not use F22s: S22 =(1/epsR) * F22s / Es , because libfreefilters.polyConj() is ony for 's' in the imaginary axis.
1173
1174
1175 ##
1176 #
1177 # This function returns the group delay as defined in [Cameron eq. 13.7, pp. 475-476] evaluated in the normalized complex frequency vector evalS .
1178 # @param Es = Characteristic polynomial of the denominator of \f$S_{21}\f$ and \f$S_{11}\f$.
1179 # @param Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$.
1180 # @param evalS = Vector with the normalized complex frequencies (s) where we want to compute the group delay.
1181 # @return groupD = Vector with the group delay at each frequency.
1182 # The group delay is computed from the formula:
1183 # \f[
1184 # \tau_{i}= Re\left[ \frac{E'(s_{i})}{E(s_{i})} - \frac{P'(s_{i})}{P(s_{i})} \right]
1185 # \f]
1186 # where \f$s_{i}\f$ are the frequency samples.
1187 #
1188 def groupDelayfunc(self, Es, Ps, evalS):
1189 groupD = np.real( Es.deriv()(evalS)/Es(evalS) - Ps.deriv()(evalS)/Ps(evalS))
1190 return groupD
1191
1192
1193 ##
1194 #
1195 # Compute S-Parameters (\f$S_{21}\f$ and \f$S_{11}\f$) from a coupling matrix.
1196 # Evaluation points are taken from attribute Sparameters.s .
1197 # Output in atributes Sparameters.S11M and Sparameters.S21M .
1198 # @param MatQ = MatrixQ class instance, containing the coupling matrix to process.
1199 # @param evalS = Normalized frequency axis (imaginary). If equal to None, it is set to 1j * np.imag(self.s). Default None.
1200 # @param checkPassive = If checkPassive is True, a warning will be raised when there are [S] parameters larger than 0 dB. Default True.
1201 # @return stErrorCM = String with a message about the error in [S] computed from M coupling matrix versus self.S11 and self.S21.
1202 # @return CM_error = Relative error in [S] computed from M coupling matrix versus self.S11 and self.S21 (float).
1203 #
1204 # By definition, \f$S_{21}(s)\f$ and \f$S_{11}(s)\f$ can be found as:
1205 # \f[
1206 # S_{21}(s) = 2\sqrt{R_{S}R_{L}} \, Z'_{N+2,1}(s)
1207 # \f]
1208 # \f[
1209 # S_{11}(s) = 1 - 2R_{S}Z'_{1,1}(s)
1210 # \f]
1211 # \f[
1212 # S_{22}(s) = 1 - 2R_{L}Z'_{N+2,N+2}(s)
1213 # \f]
1214 # where:
1215 # \f[ Z'(s)=Y'(s)^{-1} \f]
1216 # \f[ Y'(s)=jM+sI+G \f]
1217 # with \f$I \f$ a diagonal matrix with the node scaling factor squared in the diagonal except in the non-resonant nodes positions, where it is zero.
1218 # It means putting zeros in the extraNodesS first positions of the diagonal and in the last extraNodesL positions.
1219 #
1220 # \f[ G = \begin{bmatrix} G_{S} & & & & \\ & 0 & & & \\ & & \ddots & & \\ & & & 0 & \\ & & & & G_{L}\end{bmatrix} \f]
1221 # @todo Find out why sign of S11 from CM is opposite than from CP. The sign is changed below to force agreement with S11 from CP.
1222 #
1223 def fromCouplingMatrix(self, MatQ, evalS=None, checkPassive=True):
1224
1225 GS = 50/MatQ.ZoS # Normalized input load to the network
1226 GL = 50/MatQ.ZoL # Normalized output load to the network
1227
1228 # Get MatQ attributes
1229 M = MatQ.M
1230 extraNodesS = MatQ.extraNodesS
1231 extraNodesL = MatQ.extraNodesL
1232 nodeScalingFactor = MatQ.nodeScalingFactor
1233
1234 # Size of the matrix
1235 N = np.size(M,0)
1236
1237 # Evaluation points.
1238 # Since losses are now included in the diag of M, we have to evaluate in the imaginary axis
1239 if evalS is not None:
1240 eval = evalS
1241 else:
1242 eval = 1j * self.s.imag
1243
1244 # Initialize S-parameter vectors with zeros
1245 self.S11M = np.zeros(len(eval), dtype='complex128')
1246 self.S21M = np.zeros(len(eval), dtype='complex128')
1247 self.S22M = np.zeros(len(eval), dtype='complex128')
1248
1249 condY = 0 # Maximum condition number of matrices to invert
1250
1251 # For each evaluation point find S11 and S21
1252 for i in range(len(eval)):
1253 # Computation of Y
1254 I = np.diag(np.concatenate((np.zeros(extraNodesS), nodeScalingFactor[extraNodesS:-extraNodesL]**2 , np.zeros(extraNodesL))))
1255 Y = 1j*M + eval[i] * I + np.diag(np.concatenate(([GS], np.zeros(N-2) , [GL])))
1256 condY = max(condY, np.linalg.cond(Y))
1257
1258 # Computation of Z=inv(Y)
1259 Z = np.linalg.inv(Y)
1260
1261 # Computation of S11 and S22
1262 # Sign changed to agree with sign of S11 from CP
1263 self.S11M[i]= - (1.-2.*GS*Z[0, 0])
1264 self.S22M[i]= - (1.-2.*GL*Z[-1, -1])
1265
1266 # Computation of S21
1267 # Sign changed to agree with sign of S11 from CP
1268 self.S21M[i]=-2.*sqrt(GS*GL)*Z[-1, 0]
1269
1270 if condY > 10**12:
1271 warningMsg("WARNING: Computation of [S] parameters from CM:<p>Condition number of matrix to invert is %.1e.<br>[S] parameter values will be wrong due to lack of precision." % condY)
1272
1273 if checkPassive:
1274 epstol = 1e-3
1275 S11max, S21max, S22max = max(abs(self.S11M)), max(abs(self.S21M)), max(abs(self.S22M))
1276 if max([S11max, S21max, S22max]) > 1+epstol:
1277 warningMsg("""WARNING: [S] parameters computed from CM '%s' are larger than 0 dB at some frequency points<p>
1278 max[S11] = %.2f dB, max[S21] = %.2f dB, max[S22] = %.2f dB
1279 """ % (MatQ.name, 20*log10(S11max), 20*log10(S21max), 20*log10(S22max)) )
1280
1281 ##
1282 # Phase of S11 return losses computed from coupling matrix
1283 self.phaseS11M = np.unwrap(np.angle(self.S11M))
1284
1285 ##
1286 # Phase of S22 return losses computed from coupling matrix
1287 self.phaseS22M = np.unwrap(np.angle(self.S22M))
1288
1289 ##
1290 # Phase of transfer function computed from coupling matrix
1291 self.phaseS21M = np.unwrap(np.angle(self.S21M))
1292
1293 if evalS is None: # Compute group delay only for self.freq frequency axis, which has been defined when computing [S] from polynomials
1294 ##
1295 # Group delay computed from coupling matrix
1296 self.groupDelayM = -np.diff(self.phaseS21M) / (2*pi* np.diff(self.freq))
1297
1298 # Set causality principle. There may be spurious negative peaks due to computing phase at transmission zeros.
1299 self.groupDelayM[ self.groupDelayM<0 ] = 0 # CE prefers that the GD goes negative but TAS-E prefers truncation of negative values
1300
1301 ##
1302 # Frequency axis with sampling points in the middle of self.freq sampling.
1303 # It is necessary to plot self.groupDelayM
1304 self.freq2 = (self.freq[:-1] + self.freq[1:])/2
1305
1306 ##
1307 # Group delay computed from Characteristic Polynomials, sampled at self.freq2
1308 # It is necessary to compare with self.groupDelayM
1309 self.groupDelay2 = (self.groupDelay[:-1] + self.groupDelay[1:])/2
1310
1311
1312 if evalS is None: # Computation of different error measures in S11 and S21 parameters
1313
1314 errorS11dBavg = np.mean(np.abs(20.*np.log10(np.abs(self.S11)) - 20.*np.log10(np.abs(self.S11M))))
1315 errorS22dBavg = np.mean(np.abs(20.*np.log10(np.abs(self.S22)) - 20.*np.log10(np.abs(self.S22M))))
1316 errorS21dBavg = np.mean(np.abs(20.*np.log10(np.abs(self.S21)) - 20.*np.log10(np.abs(self.S21M))))
1317 errorGroupDelay = np.mean(np.abs( self.groupDelayM / self.groupDelay2 - 1 ))
1318
1319 # errorS11dBmax = np.max(np.abs(20.*np.log10(np.abs(self.S11))-20.*np.log10(np.abs(self.S11M))))
1320 # errorS11linavg = np.mean(np.abs(self.S11 - self.S11M)/np.abs(self.S11))
1321 # errorS11norm = np.linalg.norm(self.S11 - self.S11M) / np.linalg.norm(self.S11)
1322 # errorS11linmax = np.max(np.abs(self.S11 - self.S11M)/np.abs(self.S11))
1323 # myPrint("Error S11 (Avg[S11cp(dB)-S11cm(dB)]): %g dB" % errorS11dBavg )
1324 # myPrint("Error S11 (Avg[|S11cp-S11cm|/|S11cp|]): %g" % errorS11linavg )
1325 # myPrint("Error S11 ||S11cp-S11cm||/||S11cp||): %g" % errorS11norm )
1326 # myPrint("Error S11 (Max[S11cp(dB)-S11cm(dB)]): %g dB" % errorS11dBmax )
1327 # myPrint("Error S11 (Max[|S11cp-S11cm|/|S11cp|]): %g" % errorS11linmax )
1328
1329 # errorS21dBmax = np.max(np.abs(20.*np.log10(np.abs(self.S21))-20.*np.log10(np.abs(self.S21M))))
1330 # errorS21norm = np.linalg.norm(self.S21 - self.S21M) / np.linalg.norm(self.S21)
1331 # errorS21linavg = np.mean(np.abs(self.S21 - self.S21M)/np.abs(self.S21))
1332 # errorS21linmax = np.max(np.abs(self.S21 - self.S21M)/np.abs(self.S21))
1333 # myPrint("Error S21 (Avg[S21cp(dB)-S21cm(dB)]): %g dB" % errorS21dBavg )
1334 # myPrint("Error S21 (Avg[|S21cp-S11cm|/|S21cp|]): %g" % errorS21linavg )
1335 # myPrint("Error S21 ||S21cp-S21cm||/||S21cp||): %g" % errorS21norm )
1336 # myPrint("Error S21 (Max[S21cp(dB)-S21cm(dB)]): %g dB" % errorS21dBmax )
1337 # myPrint("Error S21 (Max[|S21cp-S11cm|/|S21cp|]): %g" % errorS21linmax )
1338
1339
1340 # stErrorCM = """ Average error in [S] computed from Coupling Matrix vs. [S] computed from Characteristic polynomials:
1341 # <blockquote>S<sub>11</sub> error = %.3g dB</blockquote>
1342 # <blockquote>S<sub>22</sub> error = %.3g dB</blockquote>
1343 # <blockquote>S<sub>21</sub> error = %.3g dB</blockquote>
1344 # <blockquote>Group delay mean relative error = %.9e</blockquote>""" % (errorS11dBavg, errorS22dBavg, errorS21dBavg, errorGroupDelay )
1345
1346 stErrorCM = """ Average error in [S] computed from Coupling Matrix vs. [S] computed from Characteristic polynomials:
1347 <blockquote>S<sub>11</sub> error = %.3g dB, S<sub>22</sub> error = %.3g dB, S<sub>21</sub> error = %.3g dB</blockquote>
1348 <blockquote>Group delay mean relative error = %.3e</blockquote>""" % (errorS11dBavg, errorS22dBavg, errorS21dBavg, errorGroupDelay )
1349
1350 maxErrorDB = max([errorS11dBavg, errorS21dBavg]) # Do not use errorS22dBavg, since errors in S22 may be larger
1351 CM_error = abs( 10** (-maxErrorDB/ 20) - 1 )
1352
1353 return stErrorCM, CM_error
1354 else:
1355 return None
1356
1357
1358 ##
1359 #
1360 # Compute stored energy and dissipated power at resonators and lossy couplings, from a coupling matrix.
1361 # Evaluation points are taken from attribute Sparameters.s .
1362 # Output in atributes Sparameters.energyRes, Sparameters.powerRes ,
1363 # which are 2D numpy arrays, each column is the energy(f) or power(f) at a resonator.
1364 # @param MatQ = MatrixQ class instance, containing the coupling matrix to process.
1365 #
1366 # Dissipated power at \f$i\f$-th resonator can be found as:
1367 # \f[
1368 # P_i(f) = \frac{G_i}{2} |V_i(f)|^2
1369 # \f]
1370 # being \f$ V_i \f$ the voltage at each resonator and \f$ G_i \f$ the conductance of unloaded resonators
1371 # \f[
1372 # G_i = \frac{Z_0}{Q_i \mathrm{FBW}}
1373 # \f]
1374 # where \f$ Q_i \f$ is the quality factor, \f$ \mathrm{FBW} = \Delta f / f_0 \f$ is the filter fractional bandwidth and \f$ Z_0 \f$ is the characteristic impedance of the resonators.
1375 #
1376 # The stored energy is:
1377 # \f[
1378 # W_i(f) = \frac {Q_i P_i(f)}{2 \pi f_0} = \frac{Z_0}{4 \pi f_0 \, \mathrm{FBW}} |V_i(f)|^2
1379 # \f]
1380 #
1381 # For resistive couplings between resonators \f$ i \f$ and \f$ j \f$, with conductance \f$ G_{ij} = \Im M_{ij} \f$, the voltage drop at the coupling is \f$ V_{ij}=V_i-V_j \f$ and, hence,
1382 # the dissipated power becomes:
1383 # \f[
1384 # P_{ij}(f) = \frac{G_{ij}}{2} \; |V_{ij}(f)|^2
1385 # \f]
1386 #
1387 # Following [Pozar eq. (4.14) and (4.16)], the time harmonic complex power is
1388 # \f[ P = \frac{1}{2} V I^* = \frac{1}{2} |V|^2 Y^* = Pd + 2 j \omega (W_m - W_e)
1389 # \f]
1390 # and, accounting that \f$ Y_{ij} = j M_{ij} \f$, the stored energy at magnetic or electric couplings between resonators \f$ i \f$ and \f$ j \f$,
1391 # with susceptance \f$ B_{ij} = \Im Y_{ij} = \Re M_{ij} \f$, is
1392 # \f[
1393 # W_{ij}(f) = W_m - W_e = \frac{B_{ij}}{4\omega} \; |V_{ij}(f)|^2
1394 # \f]
1395 #
1396 # Voltages at each resonator for unitary current excitation at source port are obtained as:
1397 # \f[ V(s)=Z(s) \, J(s) \f]
1398 # where \f$ J(s) \f$ is the input current vector \f$ J(s) = [I_s \;\: 0 \ldots 0]^T \f$
1399 #
1400 # If the input is matched, the input current is \f$ I_s = \sqrt{8 P_0 / Z_0} \f$. We will assume input power \f$ P_0 \f$ equal to 1 Watt.
1401 #
1402 # The impedance matrix is computed as:
1403 # \f[ Z(s)= \left( jM+sI+G \right)^{-1} \f]
1404 # with \f$I \f$ a diagonal matrix with the node scaling factor squared in the diagonal except in the non-resonant nodes positions, where it is zero.
1405 # It means putting zeros in the extraNodesS first positions of the diagonal and in the last extraNodesL positions.
1406 # \f[ G = \begin{bmatrix} G_{S} & & & & \\ & 0 & & & \\ & & \ddots & & \\ & & & 0 & \\ & & & & G_{L}\end{bmatrix} \f]
1407 #
1408 def energyPowerFromCM(self, MatQ):
1409 # Threshold to identify non-zero couplings
1410 eps = MatQ.parent.matrixElementsEps('Energy and power')
1411
1412 GS = 50/MatQ.ZoS # Normalized input load to the network
1413 GL = 50/MatQ.ZoL # Normalized output load to the network
1414
1415 # Input
1416 Po = 1. # Reference power 1W
1417 Zo = MatQ.Zo # Characteristic impedance of resonators used to compute the band-pass version of this Coupling Matrix.
1418 Is = sqrt(Po*8/Zo) # Source current
1419
1420 # Resonators Q and G
1421 Qi = MatQ.Q
1422 G = Zo / (Qi * self.FBW)
1423
1424 # Get MatQ attributes
1425 M = MatQ.M
1426 extraNodesS = MatQ.extraNodesS
1427 extraNodesL = MatQ.extraNodesL
1428 nodeScalingFactor = MatQ.nodeScalingFactor
1429
1430 # Size of the matrix
1431 N = np.size(M,0)
1432
1433 # Evaluation points.
1434 # Since losses are now included in the diag of M, we have to evaluate in the imaginary axis
1435 eval = 1j * np.imag(self.s)
1436 Nf = len(eval) # Number of frequency points
1437
1438 # Position and number of resonant nodes
1439 first = extraNodesS
1440 last = N - extraNodesL
1441 Nr = last - first
1442
1443 # Initialize output vectors with zeros
1444 energyRes, energyResRev = np.zeros((Nf, Nr)), np.zeros((Nf, Nr))
1445 powerRes, powerResRev = np.zeros((Nf, Nr)), np.zeros((Nf, Nr))
1446
1447 # Prepare data for power dissipated at resistive couplings
1448 powerCoup = []
1449 for m in range(N):
1450 for n in range(m+1, N):
1451 if n==m: continue
1452 Gmn = M[m, n].imag
1453 if np.abs(Gmn) > eps:
1454 powerCoup.append( ( (m, n), Gmn, np.zeros(Nf), np.zeros(Nf) ) )
1455
1456 # Prepare data for energy stored at magnetic or electric couplings
1457 energyCoup = []
1458 for m in range(N):
1459 for n in range(m+1, N):
1460 if n==m: continue
1461 Bmn = M[m, n].real
1462 if np.abs(Bmn) > eps:
1463 energyCoup.append( ( (m, n), Bmn, np.zeros(Nf), np.zeros(Nf) ) )
1464
1465 condY = 0 # Maximum condition number of matrices to invert
1466
1467 # For each frequency evaluation point
1468 for i in range(Nf):
1469 # Computation of Z
1470 I = np.diag(np.concatenate((np.zeros(extraNodesS), nodeScalingFactor[extraNodesS:-extraNodesL]**2 , np.zeros(extraNodesL))))
1471 Z = np.linalg.inv( 1j*M + eval[i] * I + np.diag(np.concatenate(([GS], np.zeros(N-2) , [GL]))) )
1472 condY = max(condY, np.linalg.cond(Z))
1473
1474 # Vi is the first column of Z times the current, with direct excitation, and the last column with reverse excitation
1475 V = Z[:, 0]*Is
1476 Vrev = Z[:, -1]*Is
1477
1478 V2 = np.abs(V[first:last])**2
1479 Vrev2 = np.abs(Vrev[first:last])**2
1480
1481 powerRes[i, :] = V2 * G/2
1482 powerResRev[i, :] = Vrev2 * G/2
1483
1484 energyRes[i, :] = V2 * Zo / (4*pi*self.f0*self.FBW)
1485 energyResRev[i, :] = Vrev2 * Zo / (4*pi*self.f0*self.FBW)
1486
1487 # Compute power dissipated at resistive couplings
1488 for coup in powerCoup:
1489 (m, n) = coup[0]
1490 Gmn = coup[1]
1491 coup[2][i] = np.abs(V[m] - V[n])**2 * Gmn / 2
1492 coup[3][i] = np.abs(Vrev[m] - Vrev[n])**2 * Gmn / 2
1493
1494 # Compute energy stored at magnetic or electric couplings
1495 for coup in energyCoup:
1496 (m, n) = coup[0]
1497 Bmn = coup[1]
1498 coup[2][i] = np.abs(V[m] - V[n])**2 * Bmn / (4* 2*pi*self.freq[i])
1499 coup[3][i] = np.abs(Vrev[m] - Vrev[n])**2 * Bmn / (4* 2*pi*self.freq[i])
1500
1501 if condY > 10**12:
1502 warningMsg("""WARNING: Computation of [S] parameters from CM:<p>Condition number of matrix to invert is %.1e.<br>Results will be wrong due to lack of precision.""" % condY)
1503
1504 self.powerRes = powerRes
1505 self.powerResRev = powerResRev
1506 self.energyRes = energyRes
1507 self.energyResRev = energyResRev
1508 self.powerCoup = powerCoup
1509 self.energyCoup = energyCoup
1510
1511
1512 ##
1513 #
1514 # Save [S] parameters in output file.
1515 # S11, S12, S21 and S22 are attributes of sParameters class, and of type numpy.array.
1516 # @param P = Parameter class instance.
1517 # @param prec = Number of correct significant digits to print.
1518 #
1519 def saveSparam(self, P, prec):
1520
1521 fileName = P.outDirName
1522 fout = open(fileName + '.s2p', 'w')
1523 fcsv = open(fileName + '.csv', 'w')
1524 fd = csv.writer(fcsv, delimiter = ';')
1525
1526 # Set column width
1527 colWidth = 10 + prec
1528
1529 fout.write('! %s; version %s\n' % pkgNameVersion(__doc__) )
1530 fout.write ('! [S] parameters file in TouchStone format.\n\n')
1531 fout.write('# GHZ S DB R 50\n\n')
1532
1533 fd.writerow(pkgNameVersion(__doc__))
1534 fd.writerow(['[S] parameters'])
1535
1536 st = '%-*s%-*s%-*s%-*s%-*s%-*s%-*s%-*s%-*s\n' % (colWidth, ' !Freq', colWidth, ' S11(dB)', colWidth, ' S11(ang)', colWidth, ' S12(dB)', colWidth, ' S12(ang)', colWidth, ' S21(dB)', colWidth, ' S21(ang)', colWidth, ' S22(dB)', colWidth, ' S22(ang)' )
1537 fout.write(st)
1538 fd.writerow(['Freq(GHz)', 'S11(dB)', 'S11(ang)', 'S12(dB)', 'S12(ang)', 'S21(dB)', 'S21(ang)', 'S22(dB)', 'S22(ang)'])
1539
1540 st = ''
1541 for n in range(len(self.S11)):
1542 freq = str(self.freq[n] / 1e9)
1543 S11 = complexStr(self.S11[n], prec, 'DB', colWidth)
1544 S12 = complexStr(self.S12[n], prec, 'DB', colWidth)
1545 S21 = complexStr(self.S21[n], prec, 'DB', colWidth)
1546 S22 = complexStr(self.S22[n], prec, 'DB', colWidth)
1547 st += '%-*s%-*s%-*s%-*s%-*s\n' % (colWidth, freq, colWidth, S11, colWidth, S21, colWidth, S12, colWidth, S22)
1548 fd.writerow([freq] + complexStr(self.S11[n], prec, 'DB', colWidth, listFormat = True) + complexStr(self.S12[n], prec, 'DB', colWidth, listFormat = True) + complexStr(self.S21[n], prec, 'DB', colWidth, listFormat = True) + complexStr(self.S22[n], prec, 'DB', colWidth, listFormat = True))
1549 fout.write(st)
1550 fout.close()
1551 fcsv.close()
1552
1553 myPrint("[S] parameters saved to files: %s.s2p and %s.csv" % (fileName, fileName) )
1554
1555
1556 # End class: Sparameters
1557 #!!!!!!!!!!!!!!!!!!!!!!!!!#
1558
1559 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
1560 # Begin class: M
1561
1562 ##
1563 #
1564 # Class that contains a coupling matrix, its name, the number of non-resonant nodes, the Q of resonators and some operations.
1565 #
1566 class MatrixQ(object):
1567
1568 ##
1569 #
1570 # Initialize MatrixQ class. Computes Q of resonators and stores the vector of Q values in self.Q .
1571 # @param parent = CouplingMatrices class instance to which this coupling matrix belongs.
1572 # @param M = Coupling matrix.
1573 # @param name = Matrix name.
1574 # @param extraNodesS = Number of extra rows and columns in current matrix M at the source side (non-resonant nodes).
1575 # @param extraNodesL = Number of extra rows and columns in current matrix M at the load side (non-resonant nodes).
1576 # @param FT = Frequency transformation class instance.
1577 # @param fileExt = Extension of file name, if this matrix is to be automatically saved in a file (with the output file name in the parameters class).
1578 # If equal to None, this matrix is not automatically saved. Anyway the matrix can be manually saved. Default None.
1579 # @param nodeScalingFactor = Vector containing in each node position the factor of the node scaling applied to that node. Default None.
1580 # @param flagRotateAllowed = Flag that indicates if it is allowed to do a matrix rotation (Impossible after a resonant node scaling). Default True.
1581 #
1582 def __init__(self, parent, M, name, extraNodesS, extraNodesL, FT, fileExt = None, nodeScalingFactor=None, flagRotateAllowed = True):
1583 ##
1584 # Coupling matrix
1585 self.M = M.copy() # Deep copy of numpy array
1586
1587 ##
1588 # Matrix order
1589 self.order = np.size(M, 0)
1590
1591 ##
1592 # Number of extra rows and columns in current matrix M at the source side (non-resonant nodes).
1593 self.extraNodesS = extraNodesS
1594
1595 ##
1596 # Number of extra rows and columns in current matrix M at the load side (non-resonant nodes).
1597 self.extraNodesL = extraNodesL
1598
1599 ##
1600 # Matrix name.
1601 # It is converted to unicode to allow safe comparison with matrix names entered by the user in Matrix Edit Gui.
1602 self.name = unicode(name)
1603
1604 ##
1605 # Frequency transform to use when computing Q of resonators.
1606 self.FT = FT
1607
1608 ##
1609 # Extension of file name, if this matrix is to be saved in a file. Otherwise None.
1610 self.fileExt = fileExt
1611
1612 ##
1613 # CouplingMatrices class instance to which this coupling matrix belongs.
1614 self.parent = parent
1615
1616 ##
1617 # Band pass version of the matrix
1618 self.Mbp = None
1619
1620 ##
1621 # Coupling bandwidth version of the matrix
1622 self.Mcbw = None
1623
1624 ##
1625 # Coupling bandwidth with transmission lines version of the matrix
1626 self.McbwTL = None
1627
1628 ##
1629 # Labels for rows and columns. Will be used for saving energy and power.
1630 self.rowIndices = None
1631
1632 ##
1633 # Characteristic impedance of resonators
1634 self.Zo = 50.
1635
1636 ##
1637 # Characteristic impedance of source
1638 self.ZoS = 50.
1639
1640 ##
1641 # Characteristic impedance of load
1642 self.ZoL = 50.
1643
1644 # Compute Q of resonators
1645 self.Qresonators()
1646
1647 # Vector with the information of the performed node scalings.
1648 if nodeScalingFactor is not None:
1649 self.nodeScalingFactor = nodeScalingFactor.copy()
1650 else:
1651 self.nodeScalingFactor = np.ones(np.size(self.M, 0))
1652
1653 # Flag to control if it is possible to do any matrix rotation.
1654 self.flagRotateAllowed = flagRotateAllowed
1655
1656
1657 ##
1658 #
1659 # Function to copy instances of the MatrixQ class (copies the self instance).
1660 # @return New instance of MatrixQ class copied from the self instance.
1661 #
1662 def copy(self):
1663 return MatrixQ(self.parent, self.M, self.name, self.extraNodesS, self.extraNodesL, self.FT, fileExt = self.fileExt, nodeScalingFactor= self.nodeScalingFactor, flagRotateAllowed = self.flagRotateAllowed)
1664
1665
1666 ##
1667 #
1668 # Function to set name field of the instance of the MatrixQ class.
1669 # @param name = Matrix name.
1670 #
1671 def setName(self, name):
1672 self.name = name
1673
1674
1675 ##
1676 #
1677 # Function to set file extension field of the instance of the MatrixQ class.
1678 # @param fileExt = Extension of file name, if this matrix is to be saved in a file. Otherwise None.
1679 #
1680 def setFileExt(self, fileExt):
1681 self.fileExt = fileExt
1682
1683
1684 ##
1685 #
1686 # Compute the quality factors of each resonant node of coupling matrix self.M.
1687 # Frequency transform is self.FT and output Q is stored in self.Q.
1688 #
1689 # For each node \f$ n \f$, the quality factor is
1690 # \f[
1691 # Q_{n} = \frac{1}{G_n \, \mathrm{FBW}}
1692 # \f]
1693 # where the loss conductance of the unloaded resonator is:
1694 # \f[
1695 # G_{n} = - \sum_{k} \Im M_{kn}
1696 # \f]
1697 #
1698 # <b> Proof: </b>
1699 #
1700 # In the parallel N+2 normalized lowpass prototype equivalent circuit, at each node we have current excitation \f$ I_n \f$ equal to
1701 # the current that flows away from the node through the circuit elements and couplings:
1702 # \f[
1703 # I_n = \left( G_n + j B_n + s C_n \right) V_n + \sum_{k \neq n} j M'_{nk} V_k + \sum_{k \neq n} G_{nk} (V_n - V_k)
1704 # \f]
1705 # where \f$ M'_{kn} \f$ are the real couplings, \f$ G_n \f$ is the unloaded resonator conductance (losses), \f$ G_{kn} \f$ is the conductance of resistive
1706 # couplings between resonators and, obviously, the excitation is \f$ I_n = 0 \f$ for \f$ n \neq S \f$ . For \f$ n = S \f$ we have \f$ I_n = I_S \f$ (source) .
1707 #
1708 # If we define the mutual admittances between nodes as:
1709 # \f[
1710 # Y_{nm} = \left. \frac{I_n}{V_m} \right|_{V_{k \neq m}= 0}
1711 # \f]
1712 # we have
1713 # \f[
1714 # Y_{nn} = G_n + j B_n + s C_n + \sum_{k \neq n} G_{kn}
1715 # \f]
1716 # and
1717 # \f[
1718 # Y_{nm} = j M'_{nm} - G_{nm}
1719 # \f]
1720 # In matrix form, with \f$ C_n = 1 \f$ :
1721 # \f[
1722 # [Y] = j[M] + s [I]
1723 # \f]
1724 # where the elements of the coupling matrix \f$ [M] \f$ are:
1725 # \f[
1726 # M_{nn} = \frac{Y_{nn} - s}{j} = B_n - j G_{n} - j \sum_{k \neq n} G_{kn}
1727 # \f]
1728 # and
1729 # \f[
1730 # M_{nm} = \frac{Y_{nm}}{j} = M'_{nm} + j G_{nm}
1731 # \f]
1732 # Therefore:
1733 # \f[
1734 # \Im M_{nn} = - G_{n} - \sum_{k \neq n} G_{kn} = - G_{n} - \sum_{k \neq n} \Im M_{kn}
1735 # \f]
1736 # and we can compute the resonator unloaded conductance as:
1737 # \f[
1738 # G_{n} = - \Im M_{nn} - \sum_{k \neq n} \Im M_{kn} = - \sum_{k} \Im M_{kn}
1739 # \f]
1740 #
1741 def Qresonators(self):
1742 eps = self.parent.matrixElementsEps() # Threshold to determine that losses are zero in presence of roundoff errors
1743 FBW = self.FT.FBW
1744 M = self.M
1745 first = self.extraNodesS
1746 last = np.size(M, 1) - self.extraNodesL
1747
1748 losses = np.zeros(last - first)
1749 for i in range(first, last):
1750 losses[i-first] = sum(np.imag(M[:,i]))
1751 if abs(losses[i-first]) < eps:
1752 losses[i-first] = 0.
1753
1754 ##
1755 # Q of resonators
1756 self.Q = -1 / (losses * FBW)
1757
1758
1759 ##
1760 #
1761 # Compute bandpass and coupling bandwidth versions of the matrix.
1762 #
1763 # <b>Coupling bandwidth matrix</b>
1764 #
1765 # The diagonal elements of the coupling bandwidth matrix \f$ [\Delta f] \f$ is obtained as:
1766 # \f[
1767 # \Delta f_{ii} = f_i
1768 # \f]
1769 # where the unnormalised resonant frequencies (Hz) of the resonators are:
1770 # \f[
1771 # f_i = \frac {f_0}{2} \left( \omega'' + \sqrt{ \omega'' \, ^2 + 4 } \right)
1772 # \f]
1773 # with \f$ \omega'' = - B_i \Delta f / f_0 \f$ and \f$ B_i = \Re M_{ii} \f$ . The factor \f$ \Delta f / f_0 \f$ is the fractional bandwidth.
1774 #
1775 # For non-resonant nodes \f$ f_i = f_0 \f$ .
1776 #
1777 # The couplings with source and load are:
1778 # \f[
1779 # \Delta f_{Si} = M_{Si}^2 \Delta f \qquad \Delta f_{Li} = M_{Li}^2 \Delta f
1780 # \f]
1781 # \f[
1782 # \Delta f_{SS} = \Delta f_{LL} = 0
1783 # \f]
1784 # where \f$ M_{ij} \f$ are the normalised low-pass coupling matrix elements and \f$ \Delta f \f$ is the bandwidth (MHz).
1785 #
1786 # The couplings between inner nodes (resonators or non-resonant nodes) are:
1787 # \f[
1788 # \Delta f_{ij} = M_{ij} \Delta f
1789 # \f]
1790 #
1791 # <b>Coupling bandwidth matrix with transmission lines </b>
1792 #
1793 # For each matrix element \f$ i, j \f$, the scaling factor for the pass-band coupling matrix is
1794 # \f[
1795 # K_{ij} = \alpha_{ij} M_{ij}
1796 # \f]
1797 # where \f$ [M] \f$ is the low-pass normalised matrix and \f$ [K] \f$ is the pass-band coupling matrix.
1798 # \f[
1799 # \alpha_{ij} = \frac{\pi}{2} \frac{\mathrm{FBW}}{Z_0} \qquad \mathrm{if} \; i \; \mathrm{and} \; j \neq S \; \mathrm{or} \; L
1800 # \f]
1801 # \f[
1802 # \alpha_{ij} = \sqrt{\frac{\pi}{2} \frac{\mathrm{FBW}}{Z_0 Z_{S,L}}} \qquad \mathrm{if} \; i \; \mathrm{or} \; j = S \; \mathrm{or} \; L
1803 # \f]
1804 # \f[
1805 # \alpha_{SS} = \alpha_{SL} = \alpha_{LS} = \alpha_{LL} = 1
1806 # \f]
1807 # \f[
1808 # Z_0 = Z_S = Z_L = 50 \Omega
1809 # \f]
1810 #
1811 # The diagonal elements of the coupling bandwidth matrix \f$ [\Delta f] \f$ is obtained as:
1812 # \f[
1813 # \Delta f_{ii} = f_i
1814 # \f]
1815 # where the unnormalised resonant frequencies (Hz) of the resonators are:
1816 # \f[
1817 # f_i = \frac {f_0}{2} \left( \omega'' + \sqrt{ \omega'' \, ^2 + 4 } \right)
1818 # \f]
1819 # with \f$ \omega'' = - B_i \Delta f / f_0 \f$ and \f$ B_i = \Re M_{ii} \f$ . The factor \f$ \Delta f / f_0 \f$ is the fractional bandwidth.
1820 #
1821 # For non-resonant nodes \f$ f_i = f_0 \f$ .
1822 #
1823 # The couplings with source and load are:
1824 # \f[
1825 # \Delta f_{Si} = \frac{f0}{K_{Si}} \qquad \Delta f_{Li} = \frac{f0}{K_{Li}}
1826 # \f]
1827 # \f[
1828 # \Delta f_{SS} = \Delta f_{LL} = 0
1829 # \f]
1830 #
1831 # The couplings between inner nodes (resonators or non-resonant nodes) are:
1832 # \f[
1833 # \Delta f_{ij} = f_0 \sqrt{ \left( \frac{K_{ij}}{A} \right) ^2 + C }
1834 # \f]
1835 # with
1836 # \f[
1837 # A = \frac{1}{2} \left( \frac{f_j}{f_i} + \frac{f_i}{f_j} \right)
1838 # \f]
1839 # \f[
1840 # C = \left| \frac{f_j^2 - f_i^2}{f_j^2 + f_i^2} \right| ^2
1841 # \f]
1842 #
1843 # For syncronous nodes, \f$ f_i = f_0 \f$ and therefore
1844 # \f[
1845 # A = 1 \qquad C = 0
1846 # \f]therefore
1847 # \f[
1848 # \Delta f_{ij} = K_{ij} f_0
1849 # \f]
1850 #
1851 def lp2bp2cbw(self):
1852 eps = self.parent.matrixElementsEps() # Threshold to determine that couplings are zero in presence of roundoff errors
1853
1854 N = np.size(self.M, 1)
1855 FBW = self.FT.FBW
1856 BW = self.FT.BW
1857 f0 = self.FT.f0
1858
1859 alpha = np.matrix(np.ones(N))
1860
1861 # Jordi Mateu's approach, with transmission lines
1862 alpha[0, 1:N-1] = sqrt( self.FT.FBW*pi/(2*self.Zo) ) # alpha is a 2-D matrix
1863 alpha[0, 0] = sqrt(1./self.ZoS) # ZoS is the source impedance, normally 50 ohm
1864 alpha[0, N-1] = sqrt(1./self.ZoL) # ZoL is the source impedance, normally 50 ohm
1865 self.Mbp = self.M * np.array((alpha.T*alpha))
1866
1867 # Just M*FBW, as in CMSDemo software
1868 # self.Mbp = self.M * FBW
1869
1870 self.Mcbw = np.zeros((N, N), dtype=complex)
1871 self.McbwTL = np.zeros((N, N), dtype=complex)
1872
1873 # Resonant frequencies in Hz:
1874 Bn_bp = np.diag(self.M).real * FBW
1875 fn = self.FT.unormFreq(-np.diag(self.M).real)
1876 fn2 = fn*fn
1877
1878 for m in range(N):
1879 for n in range(m+1):
1880 if abs(self.M[m, n]) < eps: continue # Null coupling
1881 if m==N-1 or n==0: # "L" row or "S" column, and their symmetrical.
1882 self.McbwTL[m, n] = (f0/1e6) / self.Mbp[m, n]
1883 self.McbwTL[n, m] = self.McbwTL[m, n]
1884 self.Mcbw[m, n] = self.M[m, n]**2 * BW / 1e6
1885 self.Mcbw[n, m] = self.Mcbw[m, n]
1886 else: # Couplings between "R" and/or "NR" nodes
1887 if n==m: # Main diagonal
1888 self.McbwTL[n, n] = fn[n] / 1e6 # MHz
1889 self.Mcbw[n, n] = fn[n] / 1e6 # MHz
1890 else: # Couplings
1891 A = ( fn[n]/fn[m] + fn[m]/fn[n] ) / 2
1892 C = abs( ( (fn2[n]-fn2[m])/(fn2[m]+fn2[m]) )**2 )
1893 # print '(m,n) =', (m, n)
1894 # print 'Kij =', self.Mbp[m, n]
1895 # print 'A, C, (Kij/A)**2, sqrt((Kij/A)**2 + C) =', A, C, (self.Mbp[m, n]/A)**2, np.sqrt( (self.Mbp[m, n]/A)**2 + abs(C) )
1896 self.McbwTL[m, n] = f0 * np.sqrt( (self.Mbp[m, n]/A)**2 + C ) / 1e6 # MHz
1897
1898 # Check if we have to take the negative sign of the sqrt, to preserve the same sign as in self.Mbp[m, n]*f0
1899 if abs(self.McbwTL[m, n] / (self.Mbp[m, n]*f0/1e6) + 1) < 1: self.McbwTL[m, n] = - self.McbwTL[m, n]
1900
1901 self.McbwTL[n, m] = self.McbwTL[m, n]
1902
1903 self.Mcbw[m, n] = self.M[m, n] * BW / 1e6
1904 self.Mcbw[n, m] = self.Mcbw[m, n]
1905
1906
1907 ##
1908 #
1909 # Rotate the coupling matrix self.M. The resulting self.M matrix maintains the eigenvalues of the original matrix and his reflection and transmission properties.
1910 # @param rotationType = Type of rotation to perform. The possible values are: {'trigonometric', 'hyperbolic'}.
1911 # @param i = Row index of the pivot that will be used in the rotation.
1912 # @param j = Column index of the pivot that will be used in the rotation.
1913 # @param rotAng = Rotation angle (rad).
1914 # @param flagQ = Flag that determines if the quality factor attribute self.Q is updated or not (True / False).
1915 #
1916 # [TN 102.2 sec. 9.1]
1917 # By definition, a rotation matrix is a matrix where every value in the diagonal is 1 except the two values \f$R_{ii}=R_{jj}=c_{r}\f$ and the rest of the values of the matrix are zero except the two values \f$R_{ji}=-R_{ij}=s_{r}\f$, where [TN 102.2 eq.(A.2) and (A.3)]:
1918 # \f[\begin{matrix}[i,j] & = & \text{rotation pivot} &\\
1919 # c_r & = & cos(rotAng), & \text{if rotationType = 'trigonometric'} \\
1920 # & & cosh(rotAng), & \text{if rotationType = 'hyperbolic'} \\
1921 # s_r & = & sin(rotAng), & \text{if rotationType = 'trigonometric'} \\
1922 # & & jsinh(rotAng), & \text{if rotationType = 'hyperbolic'} \end{matrix}
1923 # \f]
1924 #
1925 # To apply the rotation we just need to do the product [TN 102.2 eq.(A.1)]:
1926 # \f[M_{rot} = RMR^{T}\f]
1927 # where \f$R\f$ is the rotation matrix described above.
1928 #
1929 #
1930 def rotateMatrix(self, rotationType, i, j, rotAng, flagQ=True):
1931
1932 M = np.array(self.M)
1933 N = np.size(M,0) # Matrix dimension
1934
1935 if not self.flagRotateAllowed:
1936 raise synthError, u"It is not allowed to rotate matrix after a Node Scaling in a resonant node with scale factor different from 1 or -1."
1937
1938 # Row and column header labels
1939 labels = ['S']
1940 for n in range(0, self.extraNodesS-1 ): labels.append('NR')
1941 labels.extend( [ 'R' for elem in range(0, N - self.extraNodesS - self.extraNodesL) ] )
1942 for n in range(0, self.extraNodesL-1 ): labels.append('NR')
1943 labels.append('L')
1944
1945 # Check cases of forbidden rotations
1946 if i == j:
1947 raise synthError, u"It is not allowed to rotate a matrix with a pivot in the diagonal."
1948 elif labels[i] == 'S' or labels[j] == 'S' or labels[i] == 'L' or labels[j] == 'L':
1949 raise synthError, u"It is not allowed to rotate a matrix with a pivot in the S or L rows or columns."
1950 elif labels[i] != labels[j]:
1951 raise synthError, u"It is not allowed to rotate a matrix with a pivot in a row and column of different resonant/non-resonant type."
1952
1953 # Computation of the rotation angle
1954 if rotationType == 'trigonometric':
1955 cr = np.cos(rotAng)
1956 sr = np.sin(rotAng)
1957 elif rotationType == 'hyperbolic':
1958 cr = np.cosh(rotAng)
1959 sr = 1j*np.sinh(rotAng)
1960 else:
1961 raise synthError, u"Wrong rotation type. rotationType should be either 'trigonometric' or 'hyperbolic'."
1962
1963 # Construction of the rotation matrix
1964 R = np.diag(np.ones(N)+1j*np.zeros(N))
1965 R[i,i] = cr
1966 R[j,j] = cr
1967 R[i,j] = -sr
1968 R[j,i] = sr
1969
1970 # Application of the rotation matrix R to the original matrix M
1971 rotatedM = np.dot(R, np.dot(M, R.T))
1972 self.M = rotatedM
1973 if flagQ:
1974 self.Qresonators()
1975
1976
1977 # def mpRotateMatrix(self, rotationType, i, j, rotAng, flagQ=True):
1978 # """
1979 # Rotate the coupling matrix self.M, using multiprecision. The resulting self.M matrix maintains the eigenvalues of the original matrix and his reflection and transmission properties.
1980 # @param rotationType = Type of rotation to perform. The possible values are: {'trigonometric', 'hyperbolic'}.
1981 # @param i = Row index of the pivot that will be used in the rotation.
1982 # @param j = Column index of the pivot that will be used in the rotation.
1983 # @param rotAng = Rotation angle (rad).
1984 # @param flagQ = Flag that determines if the quality factor attribute self.Q is updated or not (True / False).
1985 #
1986 # [TN 102.2 sec. 9.1]
1987 # By definition, a rotation matrix is a matrix where every value in the diagonal is 1 except the two values \f$R_{ii}=R_{jj}=c_{r}\f$ and the rest of the values of the matrix are zero except the two values \f$R_{ji}=-R_{ij}=s_{r}\f$, where [TN 102.2 eq.(A.2) and (A.3)]:
1988 # \f[\begin{matrix}[i,j] & = & \text{rotation pivot} &\\
1989 # c_r & = & cos(rotAng), & \text{if rotationType = 'trigonometric'} \\
1990 # & & cosh(rotAng), & \text{if rotationType = 'hyperbolic'} \\
1991 # s_r & = & sin(rotAng), & \text{if rotationType = 'trigonometric'} \\
1992 # & & jsinh(rotAng), & \text{if rotationType = 'hyperbolic'} \end{matrix}
1993 # \f]
1994 #
1995 # To apply the rotation we just need to do the product [TN 102.2 eq.(A.1)]:
1996 # \f[M_{rot} = RMR^{T}\f]
1997 # where \f$R\f$ is the rotation matrix described above.
1998 #
1999 # """
2000 #
2001 # if not self.flagRotateAllowed:
2002 # raise synthError, u"It is not allowed to rotate matrix after a Node Scaling in a resonant node."
2003 #
2004 # M = np2mp(self.M)
2005 #
2006 # # Matrix dimension
2007 # N = np.size(self.M,0)
2008 #
2009 # # Computation of the rotation angle
2010 # if rotationType == 'trigonometric':
2011 # cr = mp.cos(rotAng)
2012 # sr = mp.sin(rotAng)
2013 # elif rotationType == 'hyperbolic':
2014 # cr = mp.cosh(rotAng)
2015 # sr = 1j*mp.sinh(rotAng)
2016 # else:
2017 # raise synthError, u"Wrong rotation type. rotationType should be either 'trigonometric' or 'hyperbolic'."
2018 #
2019 # # Construction of the rotation matrix
2020 # R = mp.eye(N)
2021 # R[i,i] = cr
2022 # R[j,j] = cr
2023 # R[i,j] = -sr
2024 # R[j,i] = sr
2025 #
2026 # # Application of the rotation matrix R to the original matrix M
2027 # rotatedM = R * M * R.T
2028 # self.M = mp2np(rotatedM)
2029 #
2030 # if flagQ:
2031 # self.Qresonators()
2032
2033
2034 ##
2035 #
2036 # Add external non-resonant nodes to the coupling matrix self.M at the source or load sides, or both.
2037 # The resulting self.M matrix will have one or two columns and rows more than the input one.
2038 # @param position = Position where we want to inflate the matrix (string). It can be 'source', 'load' or 'both'.
2039 #
2040 # If position == 'source'
2041 # \f[
2042 # Minflate=\begin{bmatrix} 0 & k_1 & 0 & \hdots \\ k_1 & & & \\ 0 & & M & \\ \vdots & & & \end{bmatrix}
2043 # \f]
2044 #
2045 # If position == 'load'
2046 # \f[
2047 # Minflate=\begin{bmatrix} & & & \vdots \\ & M & & 0 \\ & & & k_2 \\ \hdots & 0 & k_2 & 0 \end{bmatrix}
2048 # \f]
2049 #
2050 # If position == 'both'
2051 # \f[
2052 # Minflate=\begin{bmatrix} 0 & k_1 & 0 & \hdots & 0 \\ k_1 & & & & \vdots \\ 0 & & M & & 0 \\ \vdots & & & & k_2 \\ 0 & \hdots & 0 & k_2 & 0 \end{bmatrix}
2053 # \f]
2054 #
2055 # where \f$k_1\f$ and \f$k_2\f$ are usually equal to 1 unless the first or last nodes of the matrix have been previously scaled by respective values \f$k_1\f$ and \f$k_2\f$ with node-scaling operations.
2056 # If this is the case, the value of \f$k_1\f$ and \f$k_2\f$ for that nodes are retrieved from the self.nodeScalingFactor attribute.
2057 #
2058 def inflateMatrix(self, position):
2059 M = self.M
2060 # Order of the matrix
2061 N = np.size(M, 0)
2062
2063 if position == u'source':
2064 Minflate = np.zeros((N+1, N+1), dtype='complex128')
2065 Minflate[1::, 1::] = M
2066 k = self.nodeScalingFactor[0]
2067 Minflate[0, 1] = Minflate[1, 0] = k
2068 self.extraNodesS = self.extraNodesS + 1
2069 self.nodeScalingFactor = np.concatenate(([1], self.nodeScalingFactor))
2070 elif position == u'load':
2071 Minflate = np.zeros((N+1, N+1), dtype='complex128')
2072 Minflate[:-1:, :-1:] = M
2073 k = self.nodeScalingFactor[-1]
2074 Minflate[-1, -2] = Minflate[-2, -1] = k
2075 self.extraNodesL = self.extraNodesL + 1
2076 self.nodeScalingFactor = np.concatenate(( self.nodeScalingFactor, [1]))
2077 elif position == u'both':
2078 Minflate = np.zeros((N+2, N+2), dtype='complex128')
2079 Minflate[1:-1, 1:-1] = M
2080 k = self.nodeScalingFactor[0]
2081 Minflate[0, 1] = Minflate[1, 0] = k
2082 k = self.nodeScalingFactor[-1]
2083 Minflate[-1, -2] = Minflate[-2, -1] = k
2084 self.extraNodesS = self.extraNodesS + 1
2085 self.extraNodesL = self.extraNodesL + 1
2086 self.nodeScalingFactor = np.concatenate(([1], self.nodeScalingFactor, [1]))
2087 else: # Don't do anything
2088 Minflate = M
2089
2090 self.M = Minflate
2091
2092
2093 ##
2094 #
2095 # Scale a node of the coupling matrix self.M, by multiplying row [i] and column [i] by a factor. Matrix rotations are not allowed after node scaling of a resonant node.
2096 # @param i = Index of row and column which we want to scale.
2097 # @param b = Factor to use to scale the ith row and column.
2098 # @param flagQ = Flag that determines if the quality factor attribute self.Q is updated or not (True / False).
2099 #
2100 def scaleNode(self, i, b, flagQ=True):
2101 self.M[:, i] = self.M[:, i]*b
2102 self.M[i, :] = self.M[i, :]*b
2103
2104 self.nodeScalingFactor[i] = self.nodeScalingFactor[i]*b
2105 if (i >= self.extraNodesS) and (i < np.size(self.M,0) - self.extraNodesL): # If it is a resonant node
2106 if abs(b) != 1: # If the scaling factor is not equal to 1 or -1
2107 self.flagRotateAllowed = False
2108 if flagQ:
2109 self.Qresonators()
2110
2111
2112 ##
2113 #
2114 # Compute the rotation angle such that, after applying a rotation of type rotationType, the element \f$ M[m,n] \f$ of the matrix self.M is equal to zero.
2115 # The element to set to zero [m,n] must be in the row or column of the rotation pivot [i,j] or its transpose [j,i].
2116 # @param rotationType = Type of rotation to perform. The possible values are: {'trigonometric', 'hyperbolic'}.
2117 # @param i = Row index of the pivot that will be used in the rotation.
2118 # @param j = Column index of the pivot that will be used in the rotation.
2119 # @param m = Row index of the element that will become zero after rotation.
2120 # @param n = Column index of the element that will become zero after rotation.
2121 # @return rotAng = Rotation angle to remove the element \f$ M[m,n] \f$.
2122 #
2123 # First of all, and following the rules in the next table [TN 102.2 tab. A.1], we need to compute a value \f$ val \f$ which depends on the position \f$ [m,n] \f$ of the element to eliminate and
2124 # on the pivot \f$ [i,j] \f$:
2125 # \f[
2126 # \begin{tabular}{|l|l|l|l|l|l|}
2127 # \hline
2128 # & Element to eliminate & $i$ & $j$ & $k$ & $val$ \\
2129 # \hline
2130 # 1 & $M_{mn}$ & $m$ & $\neq m,n$ & $n$ & $\frac{M_{ik}}{M_{jk}}$\\
2131 # \hline
2132 # 2 & $M_{mn}$ & $\neq m,n$ & $m$ & $n$ & $-\frac{M_{jk}}{M_{ik}}$\\
2133 # \hline
2134 # 3 & $M_{mn}$ & $n$ & $\neq m,n$ & $m$ & $\frac{M_{ki}}{M_{kj}}$\\
2135 # \hline
2136 # 4 & $M_{mn}$ & $\neq m,n$ & $n$ & $m$ & $-\frac{M_{kj}}{M_{ki}}$\\
2137 # \hline
2138 # 5 & $M_{mm}$ & $m$ & $\neq m$ & - & $\frac{-M_{ij}+\sqrt{M_{ij}^{2}-M_{ii}M_{jj}}}{M_{jj}}$\\
2139 # \hline
2140 # 6 & $M_{mm}$ & $\neq m$ & $m$ & - & $\frac{M_{ij}+\sqrt{M_{ij}^{2}-M_{ii}M_{jj}}}{M_{ii}}$\\
2141 # \hline
2142 # 7 & $M_{mn}$ & $m$ & $n$ & - & $\frac{2M_{ij}}{M_{jj}-M_{ii}}$\\
2143 # \hline
2144 # 8 & $M_{mn}$ & $n$ & $m$ & - & $\frac{2M_{ij}}{M_{jj}-M_{ii}}$\\
2145 # \hline
2146 # \end{tabular}
2147 # \f]
2148 # We define \f$f_{l}=1\f$ for the cases {1,2,3,4,5,6} and \f$f_{l}=0.5\f$ for cases {7,8}.
2149 # Then, the rotation angle which sets the element \f$[m,n]\f$ to zero after applying the rotation is:
2150 # \f[\begin{matrix}
2151 # rotAng & = & f_{l}tan^{-1}(val), & \text{if rotationType = 'trigonometric'} \\
2152 # rotAng & = & -jf_{l}tan^{-1}(val), & \text{if rotationType = 'hyperbolic'} \\
2153 # \end{matrix}
2154 # \f]
2155 #
2156 def rotAngleEliminate(self, rotationType, i, j, m, n):
2157 f=1
2158 M = self.M
2159
2160 # Val computation
2161 if (i<>m) & (i<>n) & (j<>m) & (j<>n):
2162 # If the value to set to zero [m,n] is not in column i, in row i, in column j or in row j, then it cannot be removed
2163 raise synthError, u"The element to set to zero [m,n] must be in the row or column of the rotation pivot [i,j] or its transpose [j,i]."
2164 elif (m<>n) & (i==m) & (j<>m) & (j<>n):
2165 # Case 1
2166 k=n
2167 val = M[i,k]/M[j,k]
2168 elif (m<>n) & (i<>m) & (i<>n) & (j==m):
2169 # Case 2
2170 k=n
2171 val = -M[j,k]/M[i,k]
2172 elif (m<>n) & (i==n) & (j<>m) & (j<>n):
2173 # Case 3
2174 k=m
2175 val = M[k,i]/M[k,j]
2176 elif (m<>n) & (i<>m) & (i<>n) & (j==n):
2177 # Case 4
2178 k=m
2179 val = -M[k,j]/M[k,i]
2180 elif (m==n) & (i==m) & (j<>m):
2181 # Case 5
2182 if M[j,j]:
2183 val = (-M[i,j]+np.sqrt(M[i,j]**2-M[i,i]*M[j,j]))/M[j,j]
2184 else:
2185 raise synthError, u"M[j,j] is equal to zero and therefore this pivot should be changed with another one which is different from zero."
2186 elif (m==n) & (i<>m) & (j==m):
2187 # Case 6
2188 if M[i,i]:
2189 val = (M[i,j]+np.sqrt(M[i,j]**2-M[i,i]*M[j,j]))/M[i,i]
2190 else:
2191 raise synthError, u"M[i,i] is equal to zero and therefore this pivot should be changed with another one which is different from zero."
2192 else:
2193 # Cases 7 and 8
2194 f = 0.5
2195 val = 2.*M[i,j]/(M[j,j]-M[i,i])
2196
2197 # Rotation angle computation
2198 if rotationType == 'trigonometric':
2199 rotAng = f*np.arctan(val)
2200 elif rotationType == 'hyperbolic':
2201 # arctanh(val/j) = -j arctan(val)
2202 rotAng = -1j*f*np.arctan(val)
2203 else:
2204 raise synthError, u"Wrong rotation type. rotationType should be either 'trigonometric' or 'hyperbolic'."
2205
2206 return rotAng
2207
2208
2209 ##
2210 #
2211 # Add losses to resonators by substracting a factor 1/(Qeff*FBW) to the diagonal of the coupling matrix.
2212 # The typical application is to set a prescribed flatness in the ideal response of a lossy filter synthesis
2213 # by emulating the non-flat response of a lossy filter of quality factor equal to Qeff.
2214 # Qeff = Quality factor of the resonators to emulate.
2215 #
2216 def addLossesQeff(self, Qeff):
2217 Geff = 1/(Qeff*self.FT.FBW)
2218 M = self.M # Pointer to matrix, to simply code
2219 N = np.size(M, 0)
2220 first = self.extraNodesS
2221 last = N - self.extraNodesL
2222
2223 for n in range(first, last):
2224 M[n, n] = M[n, n].real + 1j * (M[n, n].imag - Geff)
2225
2226
2227 ##
2228 #
2229 # Prepare optimization parameters depending on matrix topology and forced signs.
2230 #
2231 # @param topologyReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part.
2232 # @param topologyImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part.
2233 # @param topologySignReal = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any real part, respectively.
2234 # @param topologySignImag = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any imaginary part, respectively.
2235 #
2236 # @return x0 = Initial guess. Double numpy vector with the values of the matrix entries which are different from zero.
2237 # It is a concatenation of the non-zero real elements and the non-zero imaginary elements.
2238 # @return NonZeroReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part in the upper-triangular part, including the diagonal, of the matrix.
2239 # @return NonZeroImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part in the upper-triangular part, including the diagonal, of the matrix.
2240 # The entries corresponding to resonators are False
2241 # @return numReal = Number of entries in the upper triangular part of the matrix, including the diagonal, with non-zero real part.
2242 # @return bounds = List of tuples (Min,Max) with bounds for all variables to optimize. The maximum abs of any matrix element is set to 4 times the largest coupling, with a minimum of 2.
2243 # This bound also affects the real part of diagonal elements (normalized resonant frequencies).
2244 #
2245 def setupOptimTopology(self, topologyReal, topologyImag, topologySignReal, topologySignImag):
2246
2247 M = self.M # Pointer to matrix, to simplify code
2248 N = np.size(M, 0)
2249
2250 # The matrix must be symmetric, so only the upper triangular part is considered
2251 # NonZeroReal, NonZeroImag = map(bool, np.triu(topologyReal)), map(bool, np.triu(topologyImag)) # The result of np.triu is an array of int64, we need an array of Bool
2252 UpperBool = np.fromfunction(lambda i,j: i<=j, (N, N))
2253 NonZeroReal,NonZeroImag = topologyReal*UpperBool, topologyImag*UpperBool
2254
2255 # For the resonating nodes in the diagonal, the only degree of freedom is the real part (with any sign), since the imaginary part will be fixed by the Q of resonators.
2256 Rindices = range(self.extraNodesS, N-self.extraNodesL)
2257 NonZeroReal[Rindices, Rindices] = True
2258 NonZeroImag[Rindices, Rindices] = False
2259 topologySignReal[Rindices, Rindices] = 0
2260
2261 # For the SS, LL and NR diagonal elements, the only degree of freedom is the imaginary part and must be negative
2262 NRindices = range(0, self.extraNodesS) + range(-self.extraNodesL, 0)
2263 NonZeroReal[NRindices, NRindices] = False
2264 NonZeroImag[NRindices, NRindices] = True
2265 topologySignImag[NRindices, NRindices] = -1
2266
2267 # Prepare array of forced signs for non-zero elements, to be used for optimization bounds
2268 signs = np.concatenate([ topologySignReal[NonZeroReal], topologySignImag[NonZeroImag] ])
2269
2270 # Number of non-zero real couplings, to access the 1st non-zero imag coupling after concatenation of the non-zero real and imaginary vectors
2271 numReal = len(M[NonZeroReal])
2272
2273 # Maximum abs of any matrix element is 4 times the largest coupling.
2274 # Be careful! With this implementation, this bound also affects the real part of diagonal elements (normalized resonant frequencies).
2275 maxAbsElem = 4 * np.max(np.abs(M))
2276 if maxAbsElem < 2: maxAbsElem = 2 # A reasonable minimum value is 2
2277
2278 bounds = [] # Optimization boundaries
2279 for n in range(len(signs)):
2280 lowerBound, upperBound = -maxAbsElem, maxAbsElem # Default case
2281 if signs[n] >0: lowerBound = 0 # Must be positive
2282 if signs[n] <0: upperBound = 0 # Must be negative
2283 bounds.append( (lowerBound, upperBound) ) # Imag elements
2284
2285 x0 = np.concatenate([np.real(M[NonZeroReal]), np.imag(M[NonZeroImag])]) # Optimization initial guess
2286
2287 # print 'NonZeroReal:'
2288 # print NonZeroReal
2289 # print 'topologySignReal:'
2290 # print topologySignReal
2291 # print 'NonZeroImag:'
2292 # print NonZeroImag
2293 # print 'topologySignImag:'
2294 # print topologySignImag
2295 # print 'x0, signs, Bounds =' , zip(x0, signs, bounds)
2296
2297 return x0, NonZeroReal, NonZeroImag, numReal, bounds
2298
2299
2300 ##
2301 #
2302 # Return a new MatrixQ object that is a copy of self, except for the M matrix which is created from the NonZeroValues, NonZeroReal, NonZeroImag and Q parameters.
2303 #
2304 # @param NonZeroValues = Double numpy vector with the values of the matrix entries which are different from zero.
2305 # It is a concatenation of the non-zero real elements and the non-zero imaginary elements.
2306 # @param NonZeroReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part in the upper-triangular part, including the diagonal, of the matrix.
2307 # @param NonZeroImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part in the upper-triangular part, including the diagonal, of the matrix.
2308 # The entries corresponding to resonators are False, since the imaginary part will be obtained from the prescribed Q.
2309 # @param numReal = Number of entries in the upper triangular part of the matrix, including the diagonal, with non-zero real part.
2310 # @param Q = Prescribed uniform Q for the resonators.
2311 #
2312 # @return newMatQ = New MatrixQ object
2313 #
2314 def newMatQFromNonZero(self, NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q):
2315
2316 # Create an empty matrix and fill it with the non-zero values
2317 N = self.order
2318 M = np.zeros((N, N), dtype = complex)
2319 M[NonZeroReal] = NonZeroValues[0:numReal]
2320 M[NonZeroImag] = M[NonZeroImag] + 1j*NonZeroValues[numReal::]
2321 M = M + np.triu(M, 1).T # Symmetrize the matrix
2322
2323 # We impose the Q of the resonators by fixing the imaginary part of the diagonal
2324 losses = -1/(Q*self.FT.FBW)
2325 Rindices = range(self.extraNodesS, N-self.extraNodesL) # Indices of resonant nodes
2326 for i in Rindices: M[i, i] = M[i, i] + 1j * (losses - np.sum(M[:, i].imag)) # Since now M[i,i].imag=0, then np.sum(M[:, i].imag) is equal to the sum of non-diag elements
2327
2328 return MatrixQ(self.parent, M, self.name, self.extraNodesS, self.extraNodesL, self.FT, fileExt = self.fileExt, nodeScalingFactor= self.nodeScalingFactor, flagRotateAllowed = self.flagRotateAllowed)
2329
2330
2331 ##
2332 #
2333 # Optimize the coupling matrix to achieve a prescribed uniform Q and [S] parameters as close as possible to those computed from the characteristic polynomials.
2334 # The non-zero coupling matrix entries to optimize are specified by topologyReal and topologyImag parameters,
2335 # which are set by the matrix topology widget in the GUI.
2336 #
2337 # @param P = Parameters class instance, containing filter, synthesis and sweep parameters.
2338 # @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
2339 # @param weights = Optimization weights for the difference in S11 and S21 between the computation from the coupling matrix and from polynomials:
2340 # List [ weightS11, weightS22, weightS21, weightGD ]
2341 # @param Q = Prescribed uniform Q for the resonators.
2342 # @param algor = Constrained optimization algorithm. String equal to 'L-BFG-S', 'TNC' or 'SLSQP'. See below.
2343 # @param maxIter = Maximum number of iterations, or function evaluations, depending on algorithm.
2344 # @param Nsamples = Number of samples to test. Half of them go to the pass band and the other half to the rejection band, between P.minFreq and P.maxFreq.
2345 # @param S11DR = Dynamic range of S11 and S22 to test.
2346 # @param S21DR = Dynamic range of S21 to test.
2347 # @param topologyReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part.
2348 # @param topologyImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part.
2349 # @param topologySignReal = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any real part, respectively.
2350 # @param topologySignImag = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any imaginary part, respectively.
2351 #
2352 # \section algors Constrained optimization algorithms:
2353 #
2354 # - <a href="http://en.wikipedia.org/wiki/Optimization"> http://en.wikipedia.org/wiki/Optimization</a>
2355 #
2356 # - <a href="http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html"> ttp://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html</a>
2357 #
2358 # In optimization, quasi-Newton methods (also known as variable metric methods) are algorithms for finding local maxima and minima of functions.
2359 # Quasi-Newton methods are based on Newton's method to find the stationary point of a function, where the gradient is 0.
2360 # Newton's method assumes that the function can be locally approximated as a quadratic in the region around the optimum,
2361 # and use the first and second derivatives (gradient and Hessian) to find the stationary point.
2362 # In Quasi-Newton methods the Hessian matrix of second derivatives of the function to be minimized does not need to be computed.
2363 # The Hessian is updated by analyzing successive gradient vectors instead (or approximate gradient evaluations computed using finite differences).
2364 # Quasi-Newton methods are a generalization of the secant method to find the root of the first derivative for multidimensional problems.
2365 # In multi-dimensions the secant equation is under-determined (has not a unique solution),
2366 # and quasi-Newton methods differ in how they constrain the solution, typically by adding a simple low-rank update to the current estimate of the Hessian.
2367 #
2368 # A necessary condition for optimality is that the gradient be zero.
2369 # Quasi-Newton's method need not converge unless the function has a quadratic Taylor expansion near an optimum.
2370 #
2371 # - <a href="http://en.wikipedia.org/wiki/Quasi-Newton_methods"> http://en.wikipedia.org/wiki/Quasi-Newton_methods</a>
2372 #
2373 # - <a href="http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization"> http://en.wikipedia.org/wiki/Newton%27s_method_in_optimization</a>
2374 #
2375 # \subsection lbfgsb Limited memory variation Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS-B)
2376 #
2377 # Limited memory variation of the Broyden–Fletcher–Goldfarb–Shanno (BFGS) update to approximate the inverse Hessian matrix.
2378 #
2379 # The BFGS method is one of the most popular quasi-Newton optimization methods. L-BFGS is a limited-memory version of BFGS that is particularly suited to problems with very large numbers of variables.
2380 #
2381 # References:
2382 #
2383 # - <a href="http://en.wikipedia.org/wiki/L-BFGS"> http://en.wikipedia.org/wiki/L-BFGS</a>
2384 #
2385 # - <a href="http://en.wikipedia.org/wiki/BFGS_method"> http://en.wikipedia.org/wiki/BFGS_method</a>
2386 #
2387 # - Wright, and Nocedal ‘Numerical Optimization’, 1999, pg. 198.
2388 #
2389 # - R. H. Byrd, P. Lu and J. Nocedal. A Limited Memory Algorithm for Bound Constrained Optimization, (1995), SIAM Journal on Scientific and Statistical Computing , 16, 5, pp. 1190-1208.
2390 #
2391 # - C. Zhu, R. H. Byrd and J. Nocedal. L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization (1997), ACM Transactions on Mathematical Software, Vol 23, Num. 4, pp. 550 - 560.
2392 #
2393 # \subsection slsqp Sequential Least SQuares Programming (SLSQP)
2394 #
2395 # SLSQP optimizer is a sequential least squares programming algorithm which uses the Han–Powell quasi–Newton method with a BFGS update of the B–matrix and
2396 # an L1–test function in the step–length algorithm. The optimizer uses a slightly modified version of Lawson and Hanson’s NNLS nonlinear least-squares solver.
2397 #
2398 # Sequential quadratic programming (SQP) is one of the most popular and robust algorithms for nonlinear continuous optimization.
2399 # The method is based on solving a series of subproblems designed to minimize a quadratic model of the objective subject to a linearization of the constraints.
2400 # If the problem is unconstrained, then the method reduces to Newton's method for finding a point where the gradient of the objective vanishes.
2401 # If the problem has only equality constraints, then the method is equivalent to applying Newton's method to the first-order optimality conditions, or Karush–Kuhn–Tucker conditions, of the problem.
2402 # A number of packages (including NPSOL, NLPQL, OPSYC, OPTIMA, MATLAB, and SQP) are founded on this approach.
2403 #
2404 # References:
2405 #
2406 # - <a href="http://en.wikipedia.org/wiki/Sequential_quadratic_programming"> http://en.wikipedia.org/wiki/Sequential_quadratic_programming</a>
2407 #
2408 # - Kraft D (1988) A software package for sequential quadratic programming. Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace Center — Institute for Flight Mechanics, Koln, Germany.
2409 #
2410 # \subsection tnc Bound-constrained truncated newton algorithm (TNC)
2411 #
2412 # Bound-constrained truncated newton algorithm using gradient information.
2413 #
2414 # Truncated-Newton methods are a family of methods suitable for solving large nonlinear optimization problems.
2415 # At each iteration, the current estimate of the solution is updated (i.e., a step is computed) by approximately solving the Newton equations using an iterative algorithm.
2416 # This results in a doubly iterative method: an outer iteration for the nonlinear optimization problem, and an inner iteration for the Newton equations.
2417 # The inner iteration is typically stopped or truncated before the solution to the Newton equations is obtained.
2418 # Over the past two decades, a solid convergence theory has been derived for the methods.
2419 # In addition, many algorithmic enhancements have been developed and studied, resulting in a number of publicly-available software packages.
2420 # The result has been a collection of powerful, flexible, and adaptable tools for large-scale nonlinear optimization.
2421 #
2422 # References:
2423 #
2424 # - <a href="http://js2007.free.fr/code/index.html#TNC"> http://js2007.free.fr/code/index.html#TNC</a>
2425 #
2426 # - <a href="http://iris.gmu.edu/~snash/nash/assets/TN_Survey/tn_survey.html"> http://iris.gmu.edu/~snash/nash/assets/TN_Survey/tn_survey.html</a>
2427 #
2428 # - Stephen G. Nash, “A Survey of Truncated-Newton Methods”, Journal of Computational and Applied Mathematics, 124 (2000), pp. 45-59.
2429 # <a href="http://iris.gmu.edu/~snash/papers/nash_3_99.ps"> http://iris.gmu.edu/~snash/papers/nash_3_99.ps</a>
2430 #
2431 #
2432 def uniformQ(self, P, CP, weights, Q, algor, maxIter, Nsamples, S11DR, S21DR, topologyReal, topologyImag, topologySignReal, topologySignImag):
2433
2434 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2435 # BEGIN: Helper function
2436 ##
2437 #
2438 # Return a SParameters class instance with sampling frequencies to optimize.
2439 # S11, S22, S21 and groupDelay attributes are initialized with those computed from the polynomials.
2440 #
2441 # @param P = Parameters class instance, containing filter, synthesis and sweep parameters.
2442 # @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
2443 # @param FT = FrequencyTransform class instance, containing f0 and BW parameters.
2444 # @param minFreq = Minimum frequency to optimize.
2445 # @param maxFreq = Maximum frequency to optimize.
2446 # @param Nsamples = Number of samples to test. Half of them go to the pass band and the other half to the rejection band, between minFreq and maxFreq.
2447 #
2448 # @return SPref = SParameters class instance with sampling frequencies to optimize and [S] parameters computed from CP.
2449 #
2450 def setupOptimFreq(P, CP, FT, minFreq, maxFreq, Nsamples):
2451
2452 # Frequency transform, to normalize or unnormalize frequencies
2453 f0 = FT.f0
2454 BW = FT.BW
2455
2456 # Frequency sampling points
2457 NpassBand = int(Nsamples/2)
2458 Noutband = int((Nsamples+1)/2)
2459 delta = BW/NpassBand # Passband sampling step
2460 freqInband = np.linspace(f0-BW/2+delta/2, f0+BW/2-delta/2, NpassBand) # Uniform sampling with 1st and last samples at 1/2 step (delta/2) from band transition
2461 freqOutBandLower = np.linspace(minFreq, f0-BW/2, int(Noutband/2) )
2462 freqOutBandUpper = np.linspace(f0+BW/2, maxFreq, int((Noutband+1)/2) )
2463 freq = np.concatenate([freqOutBandLower, freqInband, freqOutBandUpper])
2464 sEval = 1j * FT.normFreq(freq)
2465
2466 return Sparameters(P, CP, FT, freq)
2467 # END: Helper function
2468 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2469
2470 x0, NonZeroReal, NonZeroImag, numReal, bounds = self.setupOptimTopology(topologyReal, topologyImag, topologySignReal, topologySignImag)
2471
2472 SPref = setupOptimFreq(P, CP, self.FT, P.minFreq, P.maxFreq, Nsamples)
2473
2474 # Minimum values to test
2475 S11min = np.max( 20*np.log10(np.abs(SPref.S11)) ) - S11DR # Minimum value of S11 and S22 to test
2476 S21min = np.max( 20*np.log10(np.abs(SPref.S21)) ) - S21DR # Minimum value of S21 to test
2477
2478 # Optimization parameters
2479 args = (SPref, S11min, S21min, NonZeroReal, NonZeroImag, numReal, weights, Q) # Function arguments
2480
2481 print "\nCOUPLING MATRIX OPTIMIZATION:"
2482 if algor == 'L-BFG-S': NonZeroValues = scop.fmin_l_bfgs_b(self.uniformQOptimize, x0, fprime=None, args=args, bounds=bounds, approx_grad=True, maxfun=maxIter, iprint=1 )[0]
2483 elif algor == 'TNC': NonZeroValues = scop.fmin_tnc(self.uniformQOptimize, x0, fprime=None, args=args, bounds=bounds, approx_grad=True, maxfun=maxIter )[0]
2484 elif algor == 'SLSQP': NonZeroValues = scop.fmin_slsqp(self.uniformQOptimize, x0, fprime=None, args=args, bounds=bounds, iprint = 2, iter=maxIter)
2485 else: raise synthError, "Non valid optimization algorithm %s" % algor
2486 print ""
2487
2488 MatQOut = self.newMatQFromNonZero(NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q)
2489 self.M = MatQOut.M.copy() # Update self matrix with the optimized version
2490 self.Qresonators()
2491
2492
2493 ##
2494 #
2495 # Optimization objective function for uniformQ function. The output matrix is forced symmetric,
2496 # with positive imaginary part out of the diagonal and with uniform Q in all the resonators.
2497 # The imaginary part of the diagonal elements is determined so that the Q is constant,so it is not a variable to optimize.
2498 #
2499 # @param NonZeroValues = Double numpy vector with the values of the matrix entries which are different from zero.
2500 # It is a concatenation of the non-zero real elements and the non-zero imaginary elements.
2501 # @param SPref = SParameters class instance containing [S] computed from the characteristic polynomials. It is used as the reference to compare.
2502 # @param S11min = Minimum value of S11 and S22 to test.
2503 # @param S21min = Minimum value of S21 outside to test.
2504 # @param NonZeroReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part in the upper-triangular part, including the diagonal, of the matrix.
2505 # @param NonZeroImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part in the upper-triangular part, including the diagonal, of the matrix.
2506 # The entries corresponding to resonators are False, since the imaginary part will be obtained from the prescribed Q.
2507 # @param numReal = Number of entries in the upper triangular part of the matrix, including the diagonal, with non-zero real part.
2508 # @param weights = Optimization weights for the difference in S11(dB), S21(dB) and GD(inband) between the computation from the coupling matrix and from polynomials:
2509 # List [ weightS11, weightS22, weightS21, weightGD ]
2510 # @param Q = Prescribed uniform Q for the resonators.
2511 # @return goal = Objective function value to minimize.
2512 #
2513 def uniformQOptimize(self, NonZeroValues, SPref, S11min, S21min, NonZeroReal, NonZeroImag, numReal, weights, Q):
2514 weightS11 = weights[0]
2515 weightS22 = weights[1]
2516 weightS21 = weights[2]
2517 weightGD = weights[3]
2518
2519 # Create an empty matrix and fill it with the non-zero values
2520 MatQtest = self.newMatQFromNonZero(NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q)
2521
2522 # Computation of [S] from the testing coupling matrix
2523 SPref.fromCouplingMatrix(MatQtest, checkPassive=False)
2524 S11 = 20*np.log10(np.abs(SPref.S11))
2525 S11M = 20*np.log10(np.abs(SPref.S11M))
2526 S22 = 20*np.log10(np.abs(SPref.S22))
2527 S22M = 20*np.log10(np.abs(SPref.S22M))
2528 S21 = 20*np.log10(np.abs(SPref.S21))
2529 S21M = 20*np.log10(np.abs(SPref.S21M))
2530
2531 # Set to minimum values lower than minimum
2532 S11[S11 < S11min] = S11min
2533 S11M[S11M < S11min] = S11min
2534 S22[S22 < S11min] = S11min
2535 S22M[S22M < S11min] = S11min
2536 S21[S21 < S21min] = S21min
2537 S21M[S21M < S21min] = S21min
2538
2539 # Error between [S] from Cp and [S] from CM
2540 errS11 = np.abs(S11 - S11M)
2541 errS22 = np.abs(S22 - S22M)
2542 errS21 = np.abs( S21 - S21M)
2543 errGD = np.abs( SPref.groupDelay2 - SPref.groupDelayM ) / np.mean(np.abs(SPref.groupDelay2)) # Normalize to remove 1e-9 scale factor
2544
2545 # Objective function equaling S11 with S11M and S21 with S21M
2546 goalS11 = np.mean(errS11)
2547 goalS22 = np.mean(errS22)
2548 goalS21 = np.mean(errS21)
2549 goalGD = np.mean(errGD)
2550 goal = weightS11*goalS11 + weightS22*goalS22 + weightS21*goalS21 + weightGD*goalGD
2551
2552 return goal
2553
2554
2555 ##
2556 #
2557 # Optimize the coupling matrix to achieve a prescribed uniform Q and [S] parameters that comply with a specification mask.
2558 # The non-zero coupling matrix entries to optimize are specified by topologyReal and topologyImag parameters,
2559 # which are set by the matrix topology widget in the GUI.
2560 # Only the magnitude of [S] is optimized to comply with the S11 and S21 masks.
2561 # Group delay in the mask file, if any, is not used in optimization,
2562 #
2563 # @param P = Parameters class instance, containing filter, synthesis and sweep parameters.
2564 # @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
2565 # @param weights = Optimization weights for the difference in S11 and S21 between the computation from the coupling matrix and from polynomials:
2566 # List [ weightS11, weightS22, weightS21, weightGD ]. The last element, weightGD, is not used.
2567 # @param Q = Prescribed uniform Q for the resonators.
2568 # @param algor = Constrained optimization algorithm. String equal to 'L-BFG-S', 'TNC' or 'SLSQP'. See algorithms description MatrixQ.uniformQ method documentation.
2569 # @param maxIter = Maximum number of iterations, or function evaluations, depending on algorithm.
2570 # @param Nsamples = Number of samples to test. Half of them go to the pass band and the other half to the rejection band, between P.minFreq and P.maxFreq.
2571 # @param SM = SpecMask class instance, with the mask to compare with.
2572 # @param topologyReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part.
2573 # @param topologyImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part.
2574 # @param topologySignReal = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any real part, respectively.
2575 # @param topologySignImag = Int numpy array having -1,1 or 0 values at the position of coupling matrix entries with negative, positive or any imaginary part, respectively.
2576 #
2577 def uniformQMask(self, P, CP, weights, Q, algor, maxIter, Nsamples, SM, topologyReal, topologyImag, topologySignReal, topologySignImag):
2578
2579 x0, NonZeroReal, NonZeroImag, numReal, bounds = self.setupOptimTopology(topologyReal, topologyImag, topologySignReal, topologySignImag)
2580
2581 # Sparameters class instance, prepared to store [S] computed from CM.
2582 SPref = Sparameters(P, CP, self.FT)
2583
2584 # Mask [S] parameters to compare with
2585 SM.specToMask(SPref) # Needs SP, to get reference values of max(S21) and mean(groupDelay) for relative mask variations. The reference values are obtained from the CP.
2586 SM.maskSamplesOptimization(self.FT, Nsamples)
2587
2588 # Optimization parameters
2589 args = (SPref, SM, NonZeroReal, NonZeroImag, numReal, weights, Q) # Function arguments
2590
2591 print "\nCOUPLING MATRIX OPTIMIZATION:"
2592 if algor == 'L-BFG-S': NonZeroValues = scop.fmin_l_bfgs_b(self.uniformQOptimizeMask, x0, fprime=None, args=args, bounds=bounds, approx_grad=True, maxfun=maxIter, iprint=1 )[0]
2593 elif algor == 'TNC': NonZeroValues = scop.fmin_tnc(self.uniformQOptimizeMask, x0, fprime=None, args=args, bounds=bounds, approx_grad=True, maxfun=maxIter )[0]
2594 elif algor == 'SLSQP': NonZeroValues = scop.fmin_slsqp(self.uniformQOptimizeMask, x0, fprime=None, args=args, bounds=bounds, iprint = 2, iter=maxIter)
2595 else: raise synthError, "Non valid optimization algorithm %s" % algor
2596 print ""
2597
2598 MatQOut = self.newMatQFromNonZero(NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q)
2599 self.M = MatQOut.M.copy() # Update self matrix with the optimized version
2600 self.Qresonators()
2601
2602
2603 ##
2604 #
2605 # Optimization objective function for uniformQMask function. The output matrix is forced symmetric,
2606 # with positive imaginary part out of the diagonal and with uniform Q in all the resonators.
2607 # The imaginary part of the diagonal elements is determined so that the Q is constant,so it is not a variable to optimize.
2608 # Only the magnitude of [S] is optimized to comply with the S11 and S21 masks.
2609 # Group delay in the mask file, if any, is not used in optimization,
2610 #
2611 # @param NonZeroValues = Double numpy vector with the values of the matrix entries which are different from zero.
2612 # It is a concatenation of the non-zero real elements and the non-zero imaginary elements.
2613 # @param SPref = SParameters class instance containing method and storage to compute [S] from CM.
2614 # @param SM = SpecMask class instance, with the mask to compare with.
2615 # @param NonZeroReal = Bool numpy array having True values at the position of coupling matrix entries with non-zero real part in the upper-triangular part, including the diagonal, of the matrix.
2616 # @param NonZeroImag = Bool numpy array having True values at the position of coupling matrix entries with non-zero imaginary part in the upper-triangular part, including the diagonal, of the matrix.
2617 # The entries corresponding to resonators are False, since the imaginary part will be obtained from the prescribed Q.
2618 # @param numReal = Number of entries in the upper triangular part of the matrix, including the diagonal, with non-zero real part.
2619 # @param weights = Optimization weights for the difference in S11(dB), S21(dB) and GD(inband) between the computation from the coupling matrix and from polynomials:
2620 # List [ weightS11, weightS22, weightS21, weightGD ]. The last element, weightGD, is not used.
2621 # @param Q = Prescribed uniform Q for the resonators.
2622 # @return goal = Objective function value to minimize.
2623 #
2624 def uniformQOptimizeMask(self, NonZeroValues, SPref, SM, NonZeroReal, NonZeroImag, numReal, weights, Q):
2625
2626 weightS11 = weights[0]
2627 weightS22 = weights[1]
2628 weightS21 = weights[2]
2629
2630 # Create an empty matrix and fill it with the non-zero values
2631 MatQtest = self.newMatQFromNonZero(NonZeroValues, NonZeroReal, NonZeroImag, numReal, Q)
2632
2633 # Computation of [S] from the testing coupling matrix
2634 SPref.fromCouplingMatrix(MatQtest, evalS = SM.evalSoptim, checkPassive=False)
2635 S11M = 20*np.log10(np.abs( SPref.S11M[:SM.NfreqInBandOptim] ))
2636 S22M = 20*np.log10(np.abs( SPref.S22M[:SM.NfreqInBandOptim] ))
2637 S21MinBand = 20*np.log10(np.abs( SPref.S21M[:SM.NfreqInBandOptim] ))
2638 S21MoutBand = 20*np.log10(np.abs( SPref.S21M[SM.NfreqInBandOptim:] ))
2639
2640 # Error between [S] from Mask and [S] from CM
2641 errS11 = S11M - SM.S11optim
2642 errS11[errS11<0] = 0
2643
2644 errS22 = S22M - SM.S11optim
2645 errS22[errS22<0] = 0
2646
2647 errS21MinBand = SM.S21inBandOptim - S21MinBand
2648 errS21MinBand[errS21MinBand<0] = 0
2649
2650 errS21MoutBand = S21MoutBand - SM.S21outBandOptim
2651 errS21MoutBand[errS21MoutBand<0] = 0
2652
2653 errS21 = np.concatenate([ errS21MinBand, errS21MoutBand ])
2654
2655 # Objective function
2656 goalS11 = np.mean(errS11)
2657 goalS22 = np.mean(errS22)
2658 goalS21 = np.mean(errS21)
2659 goal = weightS11*goalS11 + weightS22*goalS22 + weightS21*goalS21
2660
2661 return goal
2662
2663
2664 ##
2665 #
2666 # Save this MatrixQ instance to an ascii file, including energy and power.
2667 # Uses self.__str__() to generate the ascii representation of the numpy matrix self.M.
2668 #
2669 # @param fileName = Name of file, including extension (string).
2670 # @param P = Parameter class instance.
2671 # @param SP = SParameters class instance.
2672 # @param precCM = Number of significant digits to save in coupling matrix and Q (int).
2673 # @param precEP = Number of significant digits to save in energy and power (int).
2674 #
2675 def saveMat(self, fileName, P, SP, precCM, precEP):
2676 #fout = codecs.open(fileName, 'w', 'utf-8') # codecs.open does not translate /n into /r/n in windows.
2677 fout = open(fileName, 'w')
2678 fout.write(u'# %s; version %s\n\n' % pkgNameVersion(__doc__) )
2679 fout.write(u'# Coupling matrix name\n')
2680 fout.write(u'Name: %s\n\n' % self.name)
2681 fout.write(u'# Number of non-resonant nodes at the source and load sides. Includes the source and load nodes.\n')
2682 fout.write(u'extraNodesS: %d\n' % self.extraNodesS)
2683 fout.write(u'extraNodesL: %d\n\n' % self.extraNodesL)
2684
2685 fout.write(u'# Vector with node scaling factors\n')
2686 st = u"nodeScalingFactor: %s\n\n" % ["%.4g" % elem for elem in self.nodeScalingFactor]
2687 fout.write(st.replace("'", ""))
2688
2689 fout.write(u'# Flag which indicates if it is allowed to do a matrix rotation (impossible after a resonant node scaling)\n')
2690 fout.write(u'flagRotateAllowed: %s\n\n' % self.flagRotateAllowed)
2691
2692 fout.write(u'# Quality factor of resonators\n')
2693 st = u"Q: %s\n\n" % ["%.4g" % elem if abs(elem) != np.inf else "%s" % str(elem) for elem in self.Q]
2694 fout.write(st.replace("'", ""))
2695
2696 fout.write(u'# Coupling matrix (each line is a matrix row):\n')
2697 st = self.__str__(precCM)
2698 fout.write(st)
2699
2700 # Compute energy and power for this matrix
2701 SP.energyPowerFromCM(self)
2702
2703 # Row and column labels
2704 size = self.M.shape[0]
2705 first = self.extraNodesS
2706 last = size - self.extraNodesL
2707 labels = ['S']
2708 for n in range(1, first ): labels.append('NR' + str(n))
2709 labels.extend( [ 'R'+str(elem) for elem in range(1, last-first+1) ] )
2710 for n in range(first, first+self.extraNodesL-1 ): labels.append('NR' + str(n))
2711 labels.append('L')
2712
2713 powerRes = SP.powerRes
2714 energyRes = SP.energyRes
2715 powerResTotal = powerRes.sum(axis=1)
2716 energyResTotal = energyRes.sum(axis=1)
2717 powerCoup =SP.powerCoup
2718 energyCoup = SP.energyCoup
2719
2720 # Set column width
2721 colWidth = 10 + precEP
2722
2723 Nf = np.size(energyRes,0) # Number of frequencies
2724 Nr = np.size(energyRes,1) # Number of resonators
2725
2726 fout.write ('\n# ENERGY AND POWER\n')
2727
2728 fout.write ('\n# Stored energy in the resonators:\n')
2729 st = '%-*s' % (colWidth, '#GHZ')
2730 for n in range(Nr): st += 'R%-*d' % (colWidth-1, n+1)
2731 st += 'Total\n'
2732 fout.write(st)
2733 st = ''
2734 for m in range(Nf):
2735 freq = str(SP.freq[m] / 1e9)
2736 st += '%-*s' % (colWidth, freq)
2737 for n in range(Nr):
2738 stNum = '%.*g' % (precEP, energyRes[m, n])
2739 st += '%-*s' % (colWidth, stNum)
2740 stNum = '%.*g' % (precEP, energyResTotal[m])
2741 st += '%-*s\n' % (colWidth, stNum)
2742 fout.write(st)
2743
2744 fout.write ('\n# Dissipated power in the resonators:\n')
2745 st = '%-*s' % (colWidth, '#GHZ')
2746 for n in range(Nr): st += 'R%-*d' % (colWidth-1, n+1)
2747 st += 'Total\n'
2748 fout.write(st)
2749 st = ''
2750 for m in range(Nf):
2751 freq = str(SP.freq[m] / 1e9)
2752 st += '%-*s' % (colWidth, freq)
2753 for n in range(Nr):
2754 stNum = '%.*g' % (precEP, powerRes[m, n])
2755 if stNum == '-0': stNum = '0'
2756 st += '%-*s' % (colWidth, stNum)
2757 stNum = '%.*g' % (precEP, powerResTotal[m])
2758 st += '%-*s\n' % (colWidth, stNum)
2759 fout.write(st)
2760
2761 if len(powerCoup) > 0:
2762 fout.write ('\n# Dissipated power in the couplings:\n')
2763 st = '%-*s' % (colWidth, '#GHZ')
2764 for coup in powerCoup:
2765 (m, n) = coup[0]
2766 elem = '%s-%s' % (labels[m], labels[n])
2767 st +='%-*s' % ( colWidth, elem)
2768 st += '\n'
2769 fout.write(st)
2770 st = ''
2771 for m in range(Nf):
2772 freq = str(SP.freq[m] / 1e9)
2773 st += '%-*s' % (colWidth, freq)
2774 for coup in powerCoup:
2775 stNum = '%.*g' % (precEP, coup[2][m])
2776 if stNum == '-0': stNum = '0'
2777 st += '%-*s' % (colWidth, stNum)
2778 st += '\n'
2779 fout.write(st)
2780
2781 fout.close()
2782
2783
2784 ##
2785 #
2786 # Convert numpy matrix to string representation.
2787 # Rounds coupling matrix elements to given precision using complexStr() function.
2788 # @param prec = Number of significant digits (int). Default 5.
2789 # @return st = String representation of the numpy matrix self.M.
2790 #
2791 def __str__(self, prec = 5):
2792
2793 M = self.M
2794 (mrows, ncols) = M.shape
2795
2796 # Set column width
2797 colWidth = 5 # Blank space between columns
2798 if (abs(M.real)>10**-prec).any(): colWidth += prec+4 # Space for the real part
2799 if (abs(M.imag)>10**-prec).any(): colWidth += prec+5 # Space for the imaginary part
2800
2801 st = '# ' + self.name + u': Low-pass matrix (normalised)\n'
2802 for m in range(mrows):
2803 for n in range(ncols):
2804 tmp = complexStr(M[m, n], prec)
2805 st += u'%-*s' % (colWidth, tmp)
2806 st += u'\n'
2807
2808 self.lp2bp2cbw()
2809
2810 st += '\n\n# ' + self.name + u': Coupling bandwidth matrix (MHz)\n'
2811 for m in range(mrows):
2812 for n in range(ncols):
2813 tmp = complexStr(self.Mcbw[m, n], prec)
2814 st += u'%-*s' % (colWidth, tmp)
2815 st += u'\n'
2816
2817 st += '\n\n# ' + self.name + u': Coupling bandwidth matrix with transmission lines (MHz)\n'
2818 for m in range(mrows):
2819 for n in range(ncols):
2820 tmp = complexStr(self.McbwTL[m, n], prec)
2821 st += u'%-*s' % (colWidth, tmp)
2822 st += u'\n'
2823
2824 return (st)
2825
2826
2827 # End class: Sparameters
2828 #!!!!!!!!!!!!!!!!!!!!!!!!!#
2829
2830
2831 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
2832 # Begin class: CouplingMatrices
2833
2834 ##
2835 #
2836 # Coupling matrices class.
2837 # If \f$ [Z] \f$ is the impedance matriz that relates the current at the network loops with the voltage sources within these loops,
2838 # the compupling matriz \f$ M \f$ is such that [TN 102.1, sec. 3.2]:
2839 # \f[
2840 # [Z] = [jM + sI]
2841 # \f]
2842 # @todo Call numpy functions as methods of array class, whenever possible.
2843 #
2844 class CouplingMatrices(object):
2845
2846 ##
2847 #
2848 # Compute Coupling Matrices for all the selected topologies.
2849 # @param P = Parameters class instance, containing filter, synthesis and sweep parameters.
2850 # @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
2851 # @param Yo = Termination admitance.
2852 # @param FT = Frequency transformation class instance.
2853 # @param SP = SParameters class instance.
2854 # @return CM = CouplingMatrices class instance.
2855 #
2856 def __init__(self, P, CP, Yo, FT, SP):
2857
2858 ##
2859 # Frequency transform used to compute these coupling matrices.
2860 # It must be used when creating new matrices or reading matrices from disk.
2861 self.FT = FT
2862
2863 ##
2864 # SParameters class instance containing the frequency sampling and SP.fromCouplingMatrix() method.
2865 # Will be used in the sensitivity analysis.
2866 self.SP = SP
2867
2868 ##
2869 # Uniform Q folded coupling matrix N+4.
2870 self.fCM_uniQ = None
2871
2872 ##
2873 # List of available coupling matrices
2874 self.listM = [ ]
2875
2876 ##
2877 # Last index of CM selection comboBox. Initially it is equal to zero(first matrix in list).
2878 self.indexCM = 0
2879
2880 ##
2881 # Roundoff error in the computation of the Tranversal Coupling Matrix N+2.
2882 # The error is obtained as the maximum of the largest of the remainders in the computation of \f$ y_{22n}(s) \f$ and \f$ y_d(s) \f$
2883 # and the constant term in the partial fraction expansion of \f$ y_{21n}(s) / Y_d(s) \f$
2884 self.CM_error = None
2885
2886 # Pre-computed coupling matrices
2887 tCM = self.transCouplingMatrix(P, CP, Yo, FT) # Transversal Coupling Matrix N+2. Definition in TN 102.1, sec. 3.2, Fig. 2b), eq. (2), (3) and (4).
2888 fCM = self.tcm2fcm(tCM) # Folded coupling matrix N+2.
2889 if CP.flagCase1: fCM_uniQ = self.uniformQFCMOrder4(fCM) # Folded coupling matrix N+4 with uniform Q.
2890 if CP.flagCase2: fCM_uniQ = self.uniformQFCMOrder6(fCM) # Folded coupling matrix N+4 with uniform Q.
2891 if CP.flagCase3: fCM_uniQ = self.uniformQFCMOrder6Case3(fCM) # Folded coupling matrix N+4 with uniform Q.
2892
2893 # Setting to zero small elements before matrix rotations can be worse because numpy sets a 0+0j where you have set a zero, and NaN appear in trigonometric operations
2894 #fCM.M [np.abs(fCM.M) < 2*self.CM_error] = 0. # Set to zero small elements of the CM, to avoid problems in matrix rotations
2895
2896 arrow = self.tcm2arrow(tCM) # arrow is never equal to None because it can always be computed.
2897 #arrow.M [np.abs(arrow.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations
2898
2899 transComplexZeros = np.array(P.transComplexZeros)
2900 transComplexZerosWithSym = np.concatenate((transComplexZeros, -transComplexZeros.conj()))
2901 normImagZeros = [ 1j*elem for elem in FT.normFreq(np.array(P.transImagZeros)) ]
2902 normComplexZeros = transComplexZerosWithSym.real + 1j * FT.normFreq(transComplexZerosWithSym.imag)
2903 TZ = np.concatenate(( normImagZeros, normComplexZeros ))
2904 if P.genChebyLP: TZ = np.concatenate([TZ, np.array(P.LPcomplexZeros).real + 1j * FT.normFreq(np.array(P.LPcomplexZeros).imag), -np.conjugate(np.array(P.LPcomplexZeros).real + 1j * FT.normFreq(np.array(P.LPcomplexZeros).imag)), np.array(P.LPrealZeros), -np.array(P.LPrealZeros) ])
2905
2906 culdesac = self.fcm2culdesac(fCM, TZ) # May be equal to None if it is not posible to compute this matrix
2907 #if culdesac:
2908 # culdesac.M [np.abs(culdesac.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations
2909
2910 triplet = self.arrow2triplet(arrow, TZ) # May be equal to None if it is not posible to compute this matrix
2911 #if triplet:
2912 # triplet.M [np.abs(triplet.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations
2913
2914 cascadedQuartets = self.fcm2cqs(fCM, 0, 'None') # May be equal to None if it is not posible to compute this matrix
2915 #if cascadedQuartets:
2916 # cascadedQuartets.M [np.abs(cascadedQuartets.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations
2917
2918 pfitzenmaier = self.fcm2pfitzenmaier(fCM) # May be equal to None if it is not posible to compute this matrix
2919 #if pfitzenmaier:
2920 # pfitzenmaier.M [np.abs(pfitzenmaier.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations
2921
2922 inlineAsym = self.fcm2inlineAsymmetric(fCM) # May be equal to None if it is not posible to compute this matrix
2923 #if inlineAsym:
2924 # inlineAsym.M [np.abs(inlineAsym.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations
2925
2926 inlineSym = self.fcm2inlineSymmetric(fCM, TZ) # May be equal to None if it is not posible to compute this matrix
2927 #if inlineSym:
2928 # inlineSym.M [np.abs(inlineSym.M) < 2*self.CM_error] = 0 # Set to zero small elements of the CM, to avoid problems in matrix rotations
2929
2930 # Append pre-computed coupling matrices to the list:
2931 self.listM.append( tCM)
2932 self.listM.append( fCM )
2933 if CP.flagCase1 or CP.flagCase2 or CP.flagCase3:
2934 self.listM.append( fCM_uniQ )
2935 self.listM.append( arrow )
2936 if culdesac: self.listM.append( culdesac )
2937 if triplet: self.listM.append( triplet )
2938 if cascadedQuartets: self.listM.append( cascadedQuartets )
2939 if pfitzenmaier: self.listM.append( pfitzenmaier )
2940 if inlineAsym: self.listM.append( inlineAsym )
2941 if inlineSym: self.listM.append( inlineSym )
2942
2943 ##
2944 # Current coupling matrix, after rotations. Is the first in the list.
2945 self.MatQ = self.listM[0]
2946
2947 ##
2948 # Backup list of self.MatQ values, for undo. Must use deep copy, since self.MatQ will be modified by edit actions.
2949 self.bkpMatQ = [ self.MatQ.copy() ]
2950
2951 ##
2952 # Backup list of actions
2953 self.bkpActions = [ ]
2954
2955
2956 ##
2957 #
2958 # Convenience function that computes the threshold to determine if coupling matrix elements (or their real or imaginary parts) are null.
2959 # @param st = String with relevant info. Default None.
2960 # @return eps = Threshold.
2961 #
2962 def matrixElementsEps(self, st = None):
2963 eps = self.CM_error*10 # Let threshold be rather close to error estimation, since elements will be filtered later in the Coupling Matrix window by using "Precision" widget.
2964 if eps < 1e-9: eps = 1e-9
2965 if eps > 1e-5: eps = 1e-5
2966
2967 if verbose and st:
2968 myPrint('Truncation threshold = %.3e for CM_error = %.3e (%s)' % (eps, self.CM_error, st) )
2969
2970 return eps
2971
2972
2973 ##
2974 #
2975 # Convenience function that updates the error estimation in self.CM_error
2976 # @param CM_error = New error estimation.
2977 # @param st = String with relevant info. Default None.
2978 #
2979 def updateCM_error(self, CM_error, st = None):
2980
2981 if verbose and st:
2982 myPrint('Updating CM error estimation %.2e with %.2e (%s)\n' % (self.CM_error, CM_error, st))
2983
2984 self.CM_error = max([self.CM_error, CM_error])
2985
2986
2987 ##
2988 #
2989 # Computation of the \f$ N+2 \f$ transversal coupling matrix of the filter defined by its caracteristic polynomials and constants.
2990 # Polynomials are of type numpy.poly1d.
2991 # This function sets the self.CM_error attribute to the roundoff error in the computation of M,
2992 # given by the largest of the remainders in the computation of \f$ y_{22n}(s) \f$ and \f$ y_d(s) \f$
2993 # and the constant term in the partial fraction expansion of \f$ y_{21n}(s) / Y_d(s) \f$
2994 # @param P = Parameters class instance, containing filter, synthesis and sweep parameters.
2995 # @param CP = CharPolys class instance, containing Characteristic Polynomials Es, Ps, Fs and constants eps, epsR.
2996 # @param Yo = Termination admitance.
2997 # @param FT = Frequency transformation class instance.
2998 # @return MatQ = MatrixQ class instance, containing the N+2 transversal coupling matrix.
2999 #
3000 # \section Smat Scattering matrix [S]
3001 #
3002 # The first step to find the \f$ N+2 \f$ transversal coupling matrix of a two-port network is finding the associated scattering matrix \f$ [S] \f$.
3003 #
3004 # By definition, the scattering matrix of a two-port network is [TN 102.1, eq (1)]:
3005 # \f[
3006 # [S]=\begin{bmatrix} S_{11}(s) & S_{12}(s) \\ S_{21}(s) & S_{22}(s) \end{bmatrix} = \frac{1}{S_{d}(s)} \begin{bmatrix} S_{11n}(s) & S_{12n}(s) \\ S_{21n}(s) & S_{22n}(s) \end{bmatrix}
3007 # \f]
3008 #
3009 # The [S] parameters are computed from the polinomials as:
3010 # \f[
3011 # S_{11} = \frac{1}{\epsilon_R} \frac{F(s)}{E(s)}
3012 # \f]
3013 # \f[
3014 # S_{21} = S_{12} = \frac{1}{\epsilon} \frac{P(s)}{E(s)}
3015 # \f]
3016 # \f[
3017 # S_{22} = \frac{1}{\epsilon_R} \frac{F_{22}(s)}{E(s)} = (-1)^N \frac{1}{\epsilon_R} \frac{F(s)^*}{E(s)}
3018 # \f]
3019 # except for asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11), where [S] parameters are computed from the lossless case above as:
3020 # \f[
3021 # S_{11} = k_{11} S_{11\mathrm{lossless}}
3022 # \f]
3023 # \f[
3024 # S_{21} = S_{12} = k_{21} S_{21\mathrm{lossless}}
3025 # \f]
3026 # \f[
3027 # S_{22} = k_{22} \left( 1 - \left( \frac{k_{21}}{k_{11}} \right) ^2 \right) + \left( \frac{k_{21}}{k_{11}} \right) ^2 S_{11 \mathrm{lossless}}
3028 # \f]
3029 #
3030 # and hence the numerators and denominator of the \f$ S_{ij} \f$ rational polynomials are:
3031 # \f[
3032 # S_{11n}(s) = \frac{F(s)}{\epsilon_{R}}
3033 # \f]
3034 # \f[
3035 # S_{12n}(s) = S_{21n} = \frac{P(s)}{\epsilon}
3036 # \f]
3037 # \f[
3038 # S_{d}(s) = E(s)
3039 # \f]
3040 # \f[
3041 # S_{22n}(s) = \frac{F_{22}(s)}{\epsilon_R}
3042 # \f]
3043 # except for asymmetrical "Prescribed Insertion Loss technique" case 1 (kS21+kS11), where they are computed as
3044 # \f[
3045 # S_{11n}(s) = k_{11} \frac{F(s)}{\epsilon_{R}}
3046 # \f]
3047 # \f[
3048 # S_{12n}(s) = S_{21n} = k_{21} \frac{P(s)}{\epsilon}
3049 # \f]
3050 # \f[
3051 # S_{d}(s) = E(s)
3052 # \f]
3053 # \f[
3054 # S_{22n}(s) = k_{22} \left( 1 - \left( \frac{k_{21}}{k_{11}} \right) ^2 \right) \, S_d(s) + \left( \frac{k_{21}}{k_{11}} \right) ^2 S_{11n}(s)
3055 # \f]
3056 # In practice, the software divides \f$ \epsilon_R \f$ by \f$ k_{11} \f$ and \f$ \epsilon \f$ by \f$ k_{21} \f$ ,
3057 # and stores the updated values to automatically apply the \f$ k_{11} \f$ and \f$ k_{21} \f$ whenever computing \f$ S_{11n} \f$ and \f$ S_{21n} \f$ .
3058 #
3059 # \section Ymat Admitance matrix [Y]
3060 #
3061 # The next step is the computation of \f$ y_{21}(s)\f$ and \f$ y_{22}(s) \f$ of the admitance matrix \f$ [Y] \f$.
3062 # [Y] is related to [S] as [TN 102.1, eq (12)]:
3063 # \f[
3064 # [Y]= \frac{1}{y_{d}(s)} \begin{bmatrix} y_{11n}(s) & y_{12n}(s) \\ y_{21n}(s) & y_{22n}(s) \end{bmatrix}
3065 # \f]
3066 # Where:
3067 # \f[
3068 # y_{21n}(s)= -2 \, S_{21n}(s)
3069 # \f]
3070 # \f[
3071 # y_{22n}(s)=\frac{(S_{d}(s)+S_{11n}(s))(S_{d}(s)-S_{22n}(s))+S_{12n}S_{21n}}{S_{d}(s)}
3072 # \f]
3073 # \f[
3074 # y_{11n}(s)=\frac{(S_{d}(s)-S_{11n}(s))(S_{d}(s)+S_{22n}(s))+S_{12n}S_{21n}}{S_{d}(s)}
3075 # \f]
3076 # \f[
3077 # y_{d}(s)=\frac{(S_{d}(s)+S_{11n}(s))(S_{d}(s)+S_{22n}(s))-S_{12n}S_{21n}}{Y_{0}(s)S_{d}(s)}
3078 # \f]
3079 # \f$ Y_{0} \f$ = Load admitance at the end of the network.
3080 #
3081 # \section Mmat N+2 Coupling Matrix
3082 # A partial fraction expansion of \f$ [y] \f$ yields [Cameron, eq. (8.41)] [TN 101.1, eq (12)]:
3083 # \f[
3084 # [Y]= \frac{1}{y_{d}(s)} \begin{bmatrix} y_{11n}(s) & y_{12n}(s) \\ y_{21n}(s) & y_{22n}(s) \end{bmatrix}
3085 # = j \begin{bmatrix} 0 & K_{\infty} \\ K_{\infty} & 0 \end{bmatrix}
3086 # + \sum_{k=1}^N \frac{1}{s-j\lambda_k} \cdot \begin{bmatrix} r_{11k}(s) & r_{12k}(s) \\ r_{21k}(s) & r_{22k}(s) \end{bmatrix}
3087 # + \begin{bmatrix} KK_{11} & 0 \\ 0 & KK_{22} \end{bmatrix}
3088 # \f]
3089 #
3090 # <ul>
3091 #
3092 # <li> Computation of the residuals \f$ r_{22k} \f$ and \f$ r_{21k} \f$:
3093 #
3094 # The residuals \f$ r_{22k} \f$ and \f$ r_{21k} \f$ can be found as [Cameron, pp. 294, footnote]:
3095 # \f[
3096 # r_{22k}(s) = \left.\frac{y_{22n}(s)}{y'_{d}(s)}\right|_{s=j\lambda_{k}} \qquad
3097 # r_{11k}(s) = \left.\frac{y_{11n}(s)}{y'_{d}(s)}\right|_{s=j\lambda_{k}} \qquad
3098 # r_{12k}(s) = r_{21k}(s) = \left.\frac{y_{21n}(s)}{y'_{d}(s)}\right|_{s=j\lambda_{k}}
3099 # \f]
3100 #
3101 # where:
3102 # <ul>
3103 # <li> \f$ y'_{d}(s) \f$ is the first derivative of \f$ y_{d}(s)\f$. It is variable Ydp in the code.
3104 # <li> \f$ j\lambda_{k} \f$ are the roots of \f$ y'_{d}(s) \f$.
3105 # </ul>
3106 # The residuals are actually computed using scipy.signal.residue() function in order to obtain the constant terms:
3107 # <ul>
3108 # <li>In \f$ [Y_{21}] \f$ the constant term is the source-load coupling term \f$ j K_{\infty} \f$ :
3109 # \f[ \frac{y_{21n}(s)}{y_d(s)} = \sum_{k=1}^N \frac{r_{21k}}{s-j\lambda_k} + j K_{\infty}
3110 # \f]
3111 # <li>In \f$ [Y_{22}] \f$ and \f$ [Y_{11}] \f$ , "Prescribed Insertion Loss technique" case 1, there is a constant term \f$ \text{\it{KK}} \not= 0 \f$
3112 # that must be included in the N+2 transversal coupling matrix [TN 102.1, eq (13)]:
3113 # \f[ \frac{y_{22n}(s)}{y_d(s)} = \sum_{k=1}^N \frac{r_{22k}}{s-j\lambda_k} + KK_{22}
3114 # \f]
3115 #
3116 # \f[ \frac{y_{11n}(s)}{y_d(s)} = \sum_{k=1}^N \frac{r_{11k}}{s-j\lambda_k} + KK_{11}
3117 # \f]
3118 # </ul>
3119 #
3120 # <li> Computation of \f$ K_{\infty} \f$:
3121 #
3122 # \f$ K_{\infty} \f$ is a real constant equal to 0 except for fully canonical filters
3123 # where the number of finite-position transmission zeros \f$ n_{fz} \f$ is equal to the filter degree \f$ N \f$.
3124 # For lossless filters it can be computed as [Cameron, eq (8.43a)]:
3125 # \f[
3126 # K_{\infty} = \frac{\epsilon_{R}}{\epsilon} \frac{1}{(\epsilon_{R}+1)}
3127 # \f]
3128 # However, this formula is not valid for lossy synthesis, so we set \f$ K_{\infty} \f$ equal to the imaginary part of the constant term in the polynomial expansion of \f$ y_{21n}'(s) / y_d(s) \f$ .
3129 #
3130 # <li> Computation of the eigenvectors \f$ T_{1k} \f$ and \f$ T_{Nk} \f$:
3131 #
3132 # The eigenvectors \f$ T_{Nk} \f$ and \f$ T_{1k} \f$ can be found as [Cameron eq. (8.55)]:
3133 # \f[
3134 # T_{Nk} = \sqrt{r_{22k}}
3135 # \f]
3136 # \f[
3137 # T_{1k} = r_{21k} / \sqrt{r_{22k}}
3138 # \f]
3139 #
3140 # <li> Build the \f$ N+2 \f$ transversal coupling matrix \f$ [M] \f$:
3141 #
3142 # The \f$ N+2 \f$ transversal coupling matrix is built from the orthogonal vectors \f$ T_{Nk} \f$ and \f$ T_{1k} \f$
3143 # and from the eigenvalues \f$ \lambda_{k} \f$, as follows [Cameron, fig 8.18]:
3144 # \f[
3145 # [M] = \begin{bmatrix}
3146 # -j KK_{11} & \cdots & T_{1k} & \cdots & K_{\infty} \\
3147 # \vdots & \ddots & & 0 & \vdots \\
3148 # T_{1k} & & -\lambda_{k} & & T_{Nk} \\
3149 # \vdots & 0 & & \ddots & \vdots \\
3150 # K_{\infty} & \cdots & T_{Nk} & \cdots & -j KK_{22}
3151 # \end{bmatrix}
3152 # \f]
3153 # The steps to build the matrix are the following:
3154 # <ol>
3155 # <li> We put the vector T1k at the first row (between the first and the last column not-included)
3156 # <li> We put the value of Kinf at the last column of the first row
3157 # <li> We put the vector TNk at the last column (between the first and the last row not-included)
3158 # <li> As the matrix M is symmetric, we copy the upper triangle to the lower
3159 # <li> We put the lambdak (with opposite sign) at the main diagonal
3160 # (between the first and the last element not-included)
3161 # <li> If there is a remaining constants \f$ KK_{11} \f$ and \f$ KK_{22} \f$
3162 # in the partial fraction expansion of \f$ Y_{11n}/Y_d \f$ and \f$ Y_{22n}/Y_d \f$,
3163 # set M[0,0] to \f$ -j KK_{11} \f$ and M[N+1,N+1] to \f$ -j KK_{22} \f$ .
3164 # </ol>
3165 # </ul>
3166 # @todo The function "r21k, p21k, remainder = sc.signal.residue(Y21np,Yd)" does not give the correct remainder when it is complex but in the tested cases it does not afect. This function should be improved for the future.
3167 #
3168 def transCouplingMatrix(self, P, CP, Yo, FT):
3169
3170 # Degree of the filtering function
3171 N = CP.Es.o
3172 # Number of finite-position transmission zeros (TZs) that have been prescribed.
3173 nfz = CP.Ps.o
3174
3175 # Scattering matrix [S] [TN 102.1, eq (1)]
3176 S11n = CP.Fs/CP.epsR
3177 S12n = CP.Ps/CP.eps
3178 S21n = CP.Ps/CP.eps
3179 Sd = CP.Es
3180
3181 if P.synthTech.title() == u'Prescribed Insertion Loss' and P.nuqTech == u'kS21+kS11':
3182 k11 = 10**(-P.nuqK11c1/20.0) # Insertion lossess
3183 k22 = 10**(-P.nuqK22c1/20.0) # Insertion lossess
3184 k21 = 10**(-P.nuqK21c1/20.0) # Insertion lossess
3185 S22n = k22 * ( 1- (k21/k11)**2) * Sd + (k21/k11)**2 * S11n
3186 else:
3187 S22n = CP.F22s/CP.epsR
3188
3189 # Admitance matrix [Y] [TN 102.1, eq (12)]:
3190 Y21n = -2*S21n
3191
3192 # Computation of Y22n and Yd with polynomial division.
3193 # This function does not work well with the example in "Prescribed_insertion_losses_K11+k21_asymmetrical.par"
3194 # (the example in TN501), while numpy.polydiv from numpy v1.3.0 works perfectly.
3195 # Y22n, error1 = self.polyExDiv( (Sd+S11n)*(Sd-S22n) + S12n*S21n, Sd)
3196 # Y11n, error2 = self.polyExDiv( (Sd-S11n)*(Sd+S22n) + S12n*S21n, Sd)
3197 # Yd, error3 = self.polyExDiv( (Sd+S11n)*(Sd+S22n) - S12n*S21n, Sd)
3198
3199 Y22n, error1 = ((Sd+S11n)*(Sd-S22n) + S12n*S21n) / Sd
3200 Y11n, error2 = ((Sd-S11n)*(Sd+S22n) + S12n*S21n) / Sd
3201 Yd, error3 = ((Sd+S11n)*(Sd+S22n) - S12n*S21n) / Sd
3202 # Although the round-off error is larger than with polynomial evaluation and polyfit, it is given by the division remainder.
3203 # numpy.polydiv and poly1d division operator fail (for numpy 1.1) with complex polynomial coeeficients.
3204
3205 if Yd.order != N: raise synthError, 'Transversal Coupling Matrix computation: Order of Yd(s) different than order of E(s).'
3206 CM_error = abs(np.concatenate( (error1.coef, error2.coef, error3.coef), 1) ).max()
3207
3208 if verbose: myPrint('\nTCM (N+2): error in polinomial division (Y22n, Y11n, Yd) = (%.2e, %.2e, %.2e)' % (abs(error1.coef).max(), abs(error2.coef).max(), abs(error3.coef).max() ))
3209
3210 # Normalization to have highest order coefficients equal to 1
3211 Y21n = Y21n/Yd.coef[0]
3212 Y22n = Y22n/Yd.coef[0]
3213 Y11n = Y11n/Yd.coef[0]
3214 Yd = Yd/Yd.coef[0]
3215
3216 # Computation of the residuals r22k and r21k [TN 102.1, eq (13)], [Cameron, pp. 294]
3217 r22k, p22k, KK22 = sc.signal.residue(Y22n,Yd)
3218 r11k, p11k, KK11 = sc.signal.residue(Y11n,Yd)
3219 #r21k, p21k, KK2 = sc.signal.residue(Y21np,Yd)
3220 r21k, p21k, KK21 = sc.signal.residue(Y21n,Yd)
3221 lambdak = p21k / 1j
3222
3223 # Kinf [Cameron, eq (8.43a)]
3224 if N<>nfz: # For non canonical filters
3225 Kinf = 0
3226 else: # For canonical filters
3227 #Kinf = CP.epsR/(CP.eps*(CP.epsR+1.)) # Only for lossless filters, does not work with lossy synthesis
3228 Kinf = KK21[-1].imag # KK22 is a polynomial (degree 0 if there are no round-off errors) and KK22[-1] is the constant term
3229
3230 # Recomputation of Y21n for canonical filters [Cameron, eq. (8.44)]. If the filter is not canonical Y21n = Y21np
3231 # It is not necessary, since we already have the correct residues.
3232 #Y21np = Y21n - 1j*Kinf*Yd
3233 #r21k, p21k, KK2 = sc.signal.residue(Y21np,Yd)
3234
3235 # Computation of the eigenvectors T1k and TNk [Cameron eq. (8.55)]:
3236 iLambdak = np.real(lambdak).argsort() # Sorting order: Decreasing real part of lambdak
3237 lambdak = lambdak[iLambdak]
3238 TNk = np.sqrt(r22k[iLambdak])
3239 T1k = r21k[iLambdak]/np.sqrt(r22k[iLambdak])
3240
3241 # Build the N+2 transversal coupling matrix [M] [Cameron, fig 8.18]:
3242 M = 1j*np.zeros((N+2,N+2))
3243
3244 # Discard constant term smaller than roundoff error
3245 if abs(KK11) <= CM_error: KK11 = np.array([0])
3246 if abs(KK22) <= CM_error: KK22 = np.array([0])
3247
3248 M[0,1:N+1] = T1k.T # 1.- We put the vector T1k at the first row (between the first and the last column not-included).
3249 M[0,N+1] = Kinf # 2.- We put the value of Kinf at the last column of the first row.
3250 M[1:N+1,N+1] = TNk # 3.- We put the vector TNk at the last column (between the first and the last row not-included)
3251 M = M + M.T # 4.- As the matrix M is symmetric, we copy the upper side to the lower
3252 # 5.- We put the lambdak (with opposite sign) at the main diagonal (between the first and the last element not-included)
3253 # with the remaining constants KK11, KK22 in the partial fraction expansion of Y22n/Yd, in M[0,0] and M[N+1,N+1]
3254 M = M - np.diag(np.concatenate((1j*KK11, lambdak, 1j*KK22), 1))
3255
3256 # If synthesis is lossless, in theory the imaginary part of couplings is roundoff error.
3257 # However, there are cases such as "Cameron - Example TCM FCM sec 8.4.3.par" in which
3258 # the residuals have imaginary parts of the order 1e-5, which are not numerical errors and are necessary
3259 # to make [S] parameters computed from CM agree [S] parameters computed from the CP.
3260
3261 if (P.synthTech.title() == u'Prescribed Insertion Loss') and (P.nuqTech == u'kS21+pole/zero'):
3262 # Remove the column and row corresponding to the (s-a)/(s+a) zero/pole that has been added in case 3 of Prescribed Insertion Loss synthesis.
3263 assert CP.nuqDelta is not None
3264
3265 # Old way to remove row and column (Txema)
3266 # eps = 1e-6
3267 # ind = np.nonzero( np.abs( - np.imag(np.diag(M)) / CP.nuqDelta - 1 ) > eps )
3268 # M = M[np.ix_(ind[0], ind[0])]
3269
3270 # New approach (JMR): Find the minimum in order to avoid problems with thresholding
3271 # We could get two poles (or no pole) closer than eps to CP.nuqDelta
3272 # It is better to get the closest pole to CP.nuqDelta
3273 indDelta = np.argmin( np.abs(np.imag(np.diag(M)) + CP.nuqDelta) )
3274 M = np.delete(np.delete(M.copy(), indDelta, 1), indDelta, 0)
3275
3276
3277 # If Minimum Insertion Loss technique with finite Q, put losses in the diagonal of M
3278 if P.synthTech.title() == u'Minimum Insertion Loss' and P.finiteQo:
3279 for n in range(N):
3280 M[n+1, n+1] += -1j / (P.Qo * FT.FBW)
3281
3282 if P.synthTech.title() == u'Predistortion':
3283 for n in range(N):
3284 M[n+1, n+1] += -1j / (P.QoPredis * FT.FBW)
3285
3286 # Set error in computation of coupling matrix
3287 self.CM_error = CM_error
3288 if verbose: myPrint('TCM (N+2): overall error = %.2e\n' % CM_error)
3289
3290 MatQ = MatrixQ(self, M, 'Transversal (N+2)', 1, 1, FT, fileExt = '.tcm')
3291
3292 # Make the sign of coupling (S,R1) equal to sign of (L,R1).
3293 # This is equivalent to changing sign of poly P(s)
3294 # This is necessary to obtain the same TCM and FCM for Cameron example at section 8.4.3 and software CMS Demo.
3295 if np.sign(MatQ.M[1,0]) != np.sign(MatQ.M[1,1]): MatQ.scaleNode(0, -1, True)
3296
3297 return MatQ
3298
3299
3300 ##
3301 #
3302 # Polynomial exact division.
3303 # Computes quot and remain polynomials so that num(s) = quot(s)*den(s) + error(s).
3304 # @param num = Dividend.
3305 # @param den = Divisor.
3306 # @return quot = Quotient.
3307 # @return error = Error in the computation of the quotient. Is equal to: num(s)-quot(s)*den(s)
3308 #
3309 # \section algPolyDiv Algorithm
3310 # Since numpy.polydiv or the equivalent numpy.poly1d polynomial division fails for complex coefficients,
3311 # at least in numpy v1.1.1, here we evaluate num(s)/den(s) for a finite number of points in the imaginary axis
3312 # and use polyfit to build the polynomial that passes trhough this points. The remainder is the error.
3313 #
3314 # This function does not work well with the example in "Prescribed_insertion_losses_K11+k21_asymmetrical.par"
3315 # (the example in journal paper publication), while numpy.polydiv from numpy v1.3.0 works perfectly.
3316 #
3317 #
3318 def polyExDiv(self, num, den):
3319 N = num.o - den.o # Order of Quotient
3320
3321 sValues = 1j * np.linspace(-2.0, 2.0, 3*(N+1)) # Use 3x the minimum number of points needed. Same result as using 40 points always.
3322 quot = np.poly1d(np.polyfit( sValues, num(sValues)/den(sValues), num.o-den.o ))
3323 error = num - quot*den
3324 return quot, error
3325
3326
3327 ##
3328 #
3329 # This function transforms the transversal \f$N+2\f$ coupling matrix \f$M\f$ into the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$.
3330 # @param MatQ = MatrixQ class instance, containing the transversal \f$N+2\f$ coupling matrix.
3331 # @return MatQfcm = MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix.
3332 #
3333 # [TN 102.2 sec. 9.2]
3334 # The procedure to obtain the matrix [FCM] consists in setting to zero step by step elements of the matrix [M] using the following:
3335 #
3336 # -Step 1: We remove all the elements of the row 1 (In total there are nerk(1) elements).
3337 #
3338 # -Step 2: We remove all the elements of the column 1 (In total there are neck(1) elements).
3339 #
3340 # -Step 3: The same as in step 1 with the row 2 (In total there are nerk(2) elements).
3341 #
3342 # -Step 4: The same as in step 2 with the column 2 (In total there are neck(2) elements).
3343 #
3344 # \f$\vdots\f$
3345 #
3346 # -Step numSteps-1: We remove all the elements of the row numRows (In total there are nerk(numRows) elements).
3347 #
3348 # -Step numSteps: We remove all the elements of the column numColumns (In total there are nerk(numColumns) elements).
3349 #
3350 # The last two steps consist in removing the last remaining row and the last remaining column of the matrix (it is not necessary that they are in that order).
3351 #
3352 def tcm2fcm(self, MatQ):
3353 Mfcm = MatQ.copy()
3354 Mfcm.setName('Folded (N+2)')
3355 Mfcm.setFileExt('.fcm')
3356 itr=0
3357 itc=0
3358
3359 # Size of the matrix
3360 N = np.size(Mfcm.M,0)
3361
3362 # Total number of steps
3363 numSteps = (N-3)*(N-2)/2
3364
3365 # Total number of rows numRows and columns numColumns to be manipulated
3366 if N%2:
3367 numRows = (N-3)/2
3368 numColumns = (N-3)/2
3369 else:
3370 numRows = (N-2)/2
3371 numColumns = (N-4)/2
3372 for k in range(numRows+numColumns):
3373 if (k+1)%2:
3374 # Elimination process of one row
3375 nerk = N-3-2*itr
3376 for l in range(nerk):
3377 m = itr
3378 n = N-l-itr-2
3379 i = N-l-itr-3
3380 j = N-l-itr-2
3381 rotAng = Mfcm.rotAngleEliminate('trigonometric', i, j, m, n)
3382 Mfcm.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3383 itr = itr+1
3384 else:
3385 # Elimination process of one column
3386 neck = N-4-2*itc
3387 for l in range(neck):
3388 m = 2+l+itc
3389 n = N-itc-1
3390 i = 2+l+itc
3391 j = 3+l+itc
3392 rotAng = Mfcm.rotAngleEliminate('trigonometric', i, j, m, n)
3393 Mfcm.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3394 itc = itc+1
3395 Mfcm.Qresonators()
3396 return Mfcm
3397
3398 ##
3399 #
3400 # This function transforms the transversal \f$N+2\f$ coupling matrix \f$M\f$ into the coupling matrix in Arrow or Wheel Form \f$Marrow\f$.
3401 # @param MatQ = MatrixQ class instance, containing the transversal \f$N+2\f$ coupling matrix.
3402 # @return MatQarrow = MatrixQ class instance, containing the Arrow topology coupling matrix.
3403 #
3404 # [Cameron]
3405 #
3406 def tcm2arrow(self, MatQ):
3407
3408 Marrow = MatQ.copy()
3409 Marrow.setName('Arrow (N+2)')
3410 Marrow.setFileExt('.arrow')
3411
3412 # Size of the matrix
3413 N = np.size(Marrow.M,0)-2
3414
3415 for i in range(1, N):
3416 for j in range(i+1, N+1):
3417 rotAng = Marrow.rotAngleEliminate('trigonometric', i, j, i-1, j)
3418 Marrow.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3419 Marrow.Qresonators()
3420 return Marrow
3421
3422 ##
3423 #
3424 # This function transforms the coupling matrix in Arrow or Wheel Form \f$Marrow\f$ into the coupling matrix in Trisections Form \f$Mtriplet\f$.
3425 # @param MatQarrow = MatrixQ class instance, containing the Arrow topology coupling matrix.
3426 # @param TZ = Transmission Zeros.
3427 # @return MatQtriplet = MatrixQ class instance, containing the Trisections topology coupling matrix.
3428 #
3429 # [Cameron]
3430 #
3431 def arrow2triplet(self, MatQarrow, TZ):
3432 Mtriplet = MatQarrow.copy()
3433 Mtriplet.setName('Triplet (N+2)')
3434 Mtriplet.setFileExt('.triplet')
3435
3436 # Size of the matrix
3437 N = np.size(Mtriplet.M,0) - 2
3438
3439 numTZ = len(TZ)
3440 for i in range(numTZ):
3441 N4 = i
3442 N3 = (N-3)-N4*2
3443 if N3>1:
3444 w01 = TZ[i]/1j
3445 rotAng = np.arctan(Mtriplet.M[N-1, N]/(w01+Mtriplet.M[N, N]))
3446 Mtriplet.rotateMatrix('trigonometric', N-1, N, rotAng, flagQ=False)
3447 for j in range(1, N3+1):
3448 rotAng = Mtriplet.rotAngleEliminate('trigonometric', N-j-1, N-j, N-j-1, N-j+1)
3449 Mtriplet.rotateMatrix('trigonometric', N-j-1, N-j, rotAng, flagQ=False)
3450 else:
3451 return None
3452
3453 Mtriplet.Qresonators()
3454 return Mtriplet
3455
3456
3457 ##
3458 #
3459 # This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Cul de Sac coupling matrix \f$Mculdesac\f$.
3460 # @param MatQfcm = MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix.
3461 # @param TZ = Transmission Zeros.
3462 # @return MatQculdesac = MatrixQ class instance, containing the Cul de Sac coupling matrix.
3463 #
3464 # [Cameron]
3465 #
3466 def fcm2culdesac(self, MatQfcm, TZ):
3467
3468 Mculdesac = MatQfcm.copy()
3469 Mculdesac.setName('Cul de Sac (N+2)')
3470 Mculdesac.setFileExt('.cds')
3471
3472 # Size of the matrix
3473 N = np.size(Mculdesac.M,0) - 2
3474
3475 numTZ = len(TZ)
3476 if N < 5 or numTZ > N-3:
3477 return None
3478
3479 if N%2:
3480 R = (N-3)/2
3481 for r in range(1, R+1):
3482 i = (N+1)/2 - r
3483 j = (N+1)/2 + r
3484 rotAng = Mculdesac.rotAngleEliminate('trigonometric', i, j, i, j-1)
3485 Mculdesac.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3486 else:
3487 R = (N-2)/2
3488 for r in range(1, R+1):
3489 i = (N+2)/2 - r
3490 j = N/2 + r
3491 rotAng = Mculdesac.rotAngleEliminate('trigonometric', i, j, i, j)
3492 Mculdesac.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3493
3494 Mculdesac.Qresonators()
3495 return Mculdesac
3496
3497 ##
3498 #
3499 # This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Cascaded Quartets coupling matrix \f$Mcqs\f$.
3500 # @param MatQfcm = MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix.
3501 # @param r = If we want to shift the quartet, r indicates the position of the first resonator node of the quartet. It should be 0 if we do not want to shift it.
3502 # @param UpDownNone = 'Up' : Shift the quartet above the diagonal. 'Down' : Shift the quartet below the diagonal. 'None' : No shifting.
3503 # @return MatQcqs = MatrixQ class instance, containing the Cascaded Quartets coupling matrix.
3504 #
3505 # [Cameron]
3506 #
3507 def fcm2cqs(self, MatQfcm, r, UpDownNone):
3508
3509 Mcqs = MatQfcm.copy()
3510 Mcqs.setName('Cascaded Quartets (N+2)')
3511 Mcqs.setFileExt('.cqs')
3512
3513 # Size of the matrix
3514 N = np.size(Mcqs.M,0) - 2
3515
3516 if N < 8:
3517 return None
3518
3519 R=(N-8)/2;
3520
3521 # Four rotations are applied
3522 i = 3+R
3523 j = 5+R
3524 A = (Mcqs.M[2+R,7+R]*Mcqs.M[3+R,4+R]*Mcqs.M[4+R,5+R]-Mcqs.M[2+R,3+R]*Mcqs.M[5+R,6+R]*Mcqs.M[6+R,7+R]+Mcqs.M[2+R,7+R]*Mcqs.M[3+R,6+R]*Mcqs.M[5+R,6+R])
3525 B = (Mcqs.M[2+R,3+R]*Mcqs.M[3+R,6+R]*Mcqs.M[6+R,7+R]-Mcqs.M[2+R,7+R]*(Mcqs.M[3+R,4+R]**2-Mcqs.M[4+R,5+R]**2-Mcqs.M[5+R,6+R]**2+Mcqs.M[3+R,6+R]**2))
3526 C = -Mcqs.M[2+R,7+R]*(Mcqs.M[3+R,6+R]*Mcqs.M[5+R,6+R]+Mcqs.M[3+R,4+R]*Mcqs.M[4+R,5+R])
3527 rotAng = np.arctan((-B+np.sqrt(B**2-4*A*C))/(2*A))
3528 Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3529
3530 i = 4+R
3531 j = 6+R
3532 rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, 3+R, j)
3533 Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3534
3535 i = 5+R
3536 j = 7+R
3537 rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, 4+R, j)
3538 Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3539
3540 i = 2+R
3541 j = 4+R
3542 rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, i, 5+R)
3543 Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3544
3545 # If we want to shift the quartet (r is the value of the shift)
3546 if UpDownNone == 'Up':
3547 # up the diagonal
3548 i = r
3549 j = r+2
3550 rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, i, r+3)
3551 Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3552 elif UpDownNone == 'Down':
3553 # down the diagonal
3554 i = r+1
3555 j = r+3
3556 rotAng = Mcqs.rotAngleEliminate('trigonometric', i, j, r, j)
3557 Mcqs.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3558 elif UpDownNone == 'None':
3559 None
3560 else:
3561 raise synthError, u"Incorrect command. UpDownNone must be 'Up', 'Down' or 'None'."
3562
3563 Mcqs.Qresonators()
3564 return Mcqs
3565
3566 ##
3567 #
3568 # This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Pfitzenmaier coupling matrix \f$Mpfi\f$.
3569 # @param MatQfcm = MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix.
3570 # @return MatQpfi = MatrixQ class instance, containing the Pfitzenmaier coupling matrix.
3571 #
3572 # [Cameron]
3573 #
3574 def fcm2pfitzenmaier(self, MatQfcm):
3575
3576 Mpfi = MatQfcm.copy()
3577 Mpfi.setName('Pfitzenmaier (N+2)')
3578 Mpfi.setFileExt('.pfi')
3579
3580 # Size of the matrix
3581 N = np.size(Mpfi.M,0) - 2
3582
3583 R = (N-4)/2
3584 for r in range(1, R+1):
3585 i = r+1
3586 j = N-i
3587 rotAng = Mpfi.rotAngleEliminate('trigonometric', i, j, i, N-r)
3588 Mpfi.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3589
3590 Mpfi.Qresonators()
3591 return Mpfi
3592
3593 ##
3594 #
3595 # This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Inline Asymmetric coupling matrix \f$Minasym\f$.
3596 # @param MatQfcm = MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix.
3597 # @return MatQinasym = MatrixQ class instance, containing the Inline Asymmetric coupling matrix.
3598 #
3599 # [Cameron]
3600 #
3601 def fcm2inlineAsymmetric(self, MatQfcm):
3602
3603 Minasym = MatQfcm.copy()
3604 Minasym.setName('Inline Asymmetric (N+2)')
3605 Minasym.setFileExt('.ina')
3606
3607 # Size of the matrix
3608 N = np.size(Minasym.M,0) - 2
3609
3610 if N == 6:
3611 i = 2
3612 j = 4
3613 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 5)
3614 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3615 elif N == 8:
3616 i = 4
3617 j = 6
3618 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 3, j)
3619 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3620 i = 2
3621 j = 4
3622 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 7)
3623 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3624 i = 3
3625 j = 5
3626 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 2, j)
3627 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3628 i = 5
3629 j = 7
3630 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 4, j)
3631 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3632 elif N == 10:
3633 i = 4
3634 j = 6
3635 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 7)
3636 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3637 i = 6
3638 j = 8
3639 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 3, j)
3640 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3641 i = 7
3642 j = 9
3643 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 6, j)
3644 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3645 elif N == 12:
3646 i = 5
3647 j = 9
3648 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 4, j)
3649 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3650 i = 3
3651 j = 5
3652 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 10)
3653 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3654 i = 2
3655 j = 4
3656 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, i, 5)
3657 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3658 i = 6
3659 j = 8
3660 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 3, j)
3661 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3662 i = 7
3663 j = 9
3664 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 6, j)
3665 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3666 i = 8
3667 j = 10
3668 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 5, j)
3669 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3670 i = 9
3671 j = 11
3672 rotAng = Minasym.rotAngleEliminate('trigonometric', i, j, 8, j)
3673 Minasym.rotateMatrix('trigonometric', i, j, rotAng, flagQ=False)
3674
3675 Minasym.Qresonators()
3676 return Minasym
3677
3678 ##
3679 #
3680 # This function transforms the coupling matrix in Folded Canonical Form (FCM) \f$Mfcm\f$ into the Inline Symmetric coupling matrix \f$Minsym\f$.
3681 # @param MatQfcm = MatrixQ class instance, containing the Folded Canonical Form (FCM) coupling matrix.
3682 # @param TZ = Transmission Zeros.
3683 # @return MatQinsym = MatrixQ class instance, containing the Inline Symmetric coupling matrix.
3684 #
3685 # [Cameron]
3686 #
3687 def fcm2inlineSymmetric(self, MatQfcm, TZ):
3688
3689 Minsym = MatQfcm.copy()
3690 Minsym.setName('Inline Symmetric (N+2)')
3691 Minsym.setFileExt('.ins')
3692
3693 # Size of the matrix
3694 N = np.size(Minsym.M,0) - 2
3695
3696 numTZ = len(TZ)
3697 if (N == 6 and numTZ == 2) or ((N == 8 or N == 10) and numTZ == 4) or (N == 12 and numTZ == 6):
3698 None
3699 else:
3700 return None
3701
3702 Meven1 = Minsym.copy()
3703 Meven1.M = Meven1.M[1:N/2+1, 1:N/2+1]
3704 Meven1.extraNodesL = Meven1.extraNodesL - 1
3705 Meven1.extraNodesS = Meven1.extraNodesS - 1
3706 Meven2 = Minsym.copy()
3707 Meven2.M = Meven2.M[1:N/2+1, N/2+1:N+1]
3708 Meven2.extraNodesL = Meven2.extraNodesL - 1
3709 Meven2.extraNodesS = Meven2.extraNodesS - 1
3710 Meven = Minsym.copy()
3711 Meven.M = Meven1.M + np.transpose(np.rot90(Meven2.M))
3712 Meven.extraNodesL = Meven.extraNodesL - 1
3713 Meven.extraNodesS = Meven.extraNodesS - 1
3714
3715 if N == 6:
3716 i = 2
3717 j = 3
3718 rotAng = np.arctan((Minsym.M[2,3]-np.sqrt(((Minsym.M[2,3])**2)-(Minsym.M[2,5]*Minsym.M[3,4])))/Minsym.M[3,4])
3719 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng, flagQ=False)
3720 elif N == 8:
3721 i = 3
3722 j = 4
3723 rotAng1 = np.arctan((Minsym.M[2,7]*Minsym.M[3,4]-np.sqrt((Minsym.M[2,7]**2)*(Minsym.M[3,4]**2)+(Minsym.M[2,7]*Minsym.M[4,5])*((Minsym.M[2,3]**2)-Minsym.M[2,7]*Minsym.M[3,6])))/(Minsym.M[2,3]**2-Minsym.M[2,7]*Minsym.M[3,6]))
3724 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng1, flagQ=False)
3725 i = 2
3726 j = 4
3727 rotAng2 = np.arctan(Minsym.M[2,7]/(Minsym.M[2,3]*np.sin(rotAng1)))
3728 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng2, flagQ=False)
3729 elif N == 10:
3730 i = 4
3731 j = 5
3732 rotAng1 = np.arctan((Minsym.M[4,5]-np.sqrt(Minsym.M[4,5]**2-Minsym.M[4,7]*Minsym.M[5,6]))/Minsym.M[5,6])
3733 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng1, flagQ=False)
3734 i = 3
3735 j = 5
3736 rotAng2 = np.arctan((np.sin(rotAng1)*Minsym.M[3,4]-np.sqrt((np.sin(rotAng1)**2*Minsym.M[3,4]**2-Minsym.M[3,8]*(Minsym.M[4,7]+Minsym.M[5,6]))))/(Minsym.M[4,7]+Minsym.M[5,6]))
3737 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng2, flagQ=False)
3738 i = 2
3739 j = 4
3740 rotAng3 = np.arctan(np.tan(rotAng2)*Minsym.M[2,3]/(Minsym.M[4,5]-np.tan(rotAng1)*Minsym.M[5,6]+np.cos(rotAng1)*np.tan(rotAng2)*Minsym.M[3,4]))
3741 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng3, flagQ=False)
3742 elif N == 12:
3743 a0 = Minsym.M[5,8]
3744 a1 = 2*Minsym.M[4,5]
3745 a2 = Minsym.M[4,9]-((Minsym.M[3,4]**2)/Minsym.M[3,10])
3746 a3 = Minsym.M[6,7]
3747 a4 = -2*Minsym.M[5,6]
3748 b0 = Minsym.M[4,9]*Minsym.M[6,7]
3749 b1 = -2*Minsym.M[4,5]*Minsym.M[6,7]
3750 b2 = Minsym.M[5,8]*Minsym.M[6,7]-Minsym.M[5,6]**2
3751 b3 = -2*Minsym.M[5,6]*Minsym.M[4,5]
3752 b4 = Minsym.M[4,9]*Minsym.M[5,8]-Minsym.M[4,5]**2
3753 b5 = 2*Minsym.M[5,6]*Minsym.M[4,9]
3754 c=a3*b5-a4*b4
3755 c0 = (a0*b4-a3*b0)/c
3756 c1 = (a1*b4-a3*b1)/c
3757 c2 = (a2*b4-a3*b2)/c
3758 c3 = a3*b3/c
3759 d = a2*c3**2+a3*c2**2
3760 d0 = (a0+a3*c0**2+a4*c0)/d
3761 d1 = (a1+2*a0*c3+2*a3*c0*c1+a4*(c1+c0*c3))/d
3762 d2 = (a2+2*a1*c3+a0*c3**2+a3*(c1**2+2*c0*c2)+a4*(c2+c1*c3))/d
3763 d3 = (2*a2*c3+a1*c3**2+2*a3*c1*c2+a4*c2*c3)/d
3764 t1 = np.roots([1., d3, d2, d1, d0])[0]
3765 i = 4
3766 j = 5
3767 rotAng1 = np.arctan(t1)
3768 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng1, flagQ=False)
3769 i = 5
3770 j = 6
3771 rotAng2 = np.arctan((c0+c1*t1+c2*t1**2)/((1+c3*t1)*np.sqrt(1+t1**2)))
3772 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng2, flagQ=False)
3773 i = 4
3774 j = 6
3775 rotAng3 = np.arctan(Meven.M[i-1, j-1]/Meven.M[j-1, j-1])
3776 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng3, flagQ=False)
3777 i = 3
3778 j = 5
3779 rotAng4 = np.arctan(Meven.M[i-1, j-1]/Meven.M[j-1, j-1])
3780 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng4, flagQ=False)
3781 i = 2
3782 j = 4
3783 rotAng5 = np.arctan(Meven.M[i-1, j]/Meven.M[j-1, j])
3784 Meven.rotateMatrix('trigonometric', i-1, j-1, rotAng5, flagQ=False)
3785
3786 Minline = np.zeros([N/2,N], dtype='complex128')
3787 # Upper diagonal
3788 for i in range(N/2):
3789 if i<N/2-1:
3790 Minline[i,i+1] = Meven.M[i,i+1]
3791 Minline[i+1,i] = Minline[i,i+1]
3792 else:
3793 Minline[i,i+1] = Meven.M[i,i]
3794
3795 # Other elements
3796 if N == 6:
3797 Minline[0,3] = Meven.M[0,2]
3798 Minline[2,5] = Minline[0,3]
3799 elif N == 8:
3800 Minline[0,3] = Meven.M[0,3]
3801 Minline[3,0] = Minline[0,3]
3802 Minline[2,5] = Meven.M[2,2]
3803 elif N == 10:
3804 Minline[0,3] = Meven.M[0,3]
3805 Minline[3,0] = Minline[0,3]
3806 Minline[2,5] = Meven.M[2,4]
3807 Minline[4,7] = Minline[2,5]
3808 elif N == 12:
3809 Minline[0,3] = Meven.M[0,3]
3810 Minline[3,0] = Minline[0,3]
3811 Minline[2,5] = Meven.M[2,5]
3812 Minline[4,7] = Meven.M[4,4]
3813 Minline[5,2] = Minline[2,5]
3814
3815 # Fulfill the two lower quadrants of the matrix, which are symmetric to the upper ones
3816 Minline1=np.rot90(Minline)
3817 Minline2=np.rot90(Minline1)
3818 Minsym.M[1:N+1, 1:N+1]=np.concatenate([Minline, Minline2])
3819
3820 Minsym.Qresonators()
3821 return Minsym
3822
3823 ##
3824 #
3825 # Read coupling matrix from disk file.
3826 # @param fileName = File name (string).
3827 # @return MatQ = Coupling matrix (MatrixQ class instance.)
3828 #
3829 def readCouplingMatrix(self, fileName):
3830 try:
3831 file = codecs.open(fileName, 'r', 'utf-8')
3832 except IOError, err:
3833 raise parseError, 'I/O error: %s' % err
3834 stList = file.readlines()
3835 file.close()
3836
3837 validLines = [ st for st in stList if not st.startswith('#') and not st.isspace() ]
3838
3839 pat = re.compile('^Name:(.+)')
3840 name = pat.match(validLines[0]).groups(1)[0].lstrip().rstrip()
3841
3842 pat = re.compile('^extraNodesS:(.+)')
3843 extraNodesS = int(pat.match(validLines[1]).groups(1)[0])
3844
3845 pat = re.compile('^extraNodesL:(.+)')
3846 extraNodesL = int(pat.match(validLines[2]).groups(1)[0])
3847
3848 pat = re.compile('^nodeScalingFactor:(.+)')
3849 nodeScalingFactor = np.array( eval(pat.match(validLines[3]).groups(1)[0].strip()) )
3850
3851 pat = re.compile('^flagRotateAllowed:(.+)')
3852 flagRotateAllowed = eval(pat.match(validLines[4]).groups(1)[0].strip())
3853
3854 # Matrix order
3855 N = len(nodeScalingFactor)
3856
3857 Mlist = []
3858 for line in validLines[6:6+N]:
3859 Mlist.append([complex(elem) for elem in line.split()])
3860 M = np.array(Mlist)
3861
3862 MatQ = MatrixQ(self, M, name, extraNodesS, extraNodesL, self.FT, nodeScalingFactor=nodeScalingFactor, flagRotateAllowed=flagRotateAllowed)
3863 myPrint("Coupling matrix '%s' read from file: %s" % (name, os.path.basename(fileName)) )
3864 return MatQ
3865
3866
3867 ##
3868 #
3869 # Save coupling matrices in output file.
3870 # The matrices are attributes of CouplingMatrices class, and of type numpy.array.
3871 # @param P = Parameter class instance.
3872 # @param SP = SParameters class instance.
3873 # @param precCM = Number of significant digits to save in coupling matrix and Q (int).
3874 # @param precEP = Number of significant digits to save in energy and power (int).
3875 #
3876 def saveCouplingMatrices(self, P, SP, precCM, precEP):
3877 # Transversal coupling matrix
3878 myPrint('Coupling matrices:')
3879 for MatQ in self.listM:
3880 if MatQ.fileExt is None: continue
3881 myPrint(MatQ.name + ':')
3882 st = " Q of resonators = %s" % ["%.4g" % elem if abs(elem) != np.inf else "%s" % str(elem) for elem in MatQ.Q] # ["%.4g" % elem for elem in MatQ.Q]
3883 myPrint(st.replace("'", ""))
3884 fileName = P.outDirName + MatQ.fileExt
3885 MatQ.saveMat(fileName, P, SP, precCM, precEP)
3886 myPrint(" Matrix saved to file: %s" % fileName)
3887
3888
3889 ##
3890 #
3891 # This function transforms the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix of a filter of order 4 into an N+4 matrix so that all the resonant nodes have equal quality factor Q.
3892 # (Valid for Butterworth, Chevyshev, Quasieliptic and Generalized Chebyshev with symmetric transmission zeros).
3893 # @param MatQfcm2 = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix.
3894 # @return MatQfcm4 = MatrixQ class instance, containing the Folded \f$N+4\f$ Canonical Form (FCM) coupling matrix with uniform Q.
3895 #
3896 # For a better description see TN102.3 sec. 3.2.
3897 #
3898 # First of all we add two non-resonant nodes, what means converting the N+2 folded matrix into the N+4. Then we apply two hyperbolic rotations so that we have a uniform Q. For doing so, it can be proved that the necessary angle is (TN102.3 eq. (10)):
3899 # \f[
3900 # \alpha = \frac{1}{2}ln\left[\frac{\frac{M_{01}^{2}}{G}\pm \sqrt{4G_{1}\frac{M_{01}^{2}}{G}-4G_{1}^{2}+16M_{12}^{2}}}{2G_{1}+\frac{M_{01}^{2}}{G}+4M_{12}}\right]
3901 # \f]
3902 # Finally, we need to re-scale with a value \f$h\f$ so that the non-resonant nodes have a sum of imaginary parts equal to 0.
3903 #
3904 def uniformQFCMOrder4(self, MatQfcm2):
3905 MatQ = MatQfcm2.copy()
3906 MatQ.setName('Folded with uniform Q (N+4)')
3907 MatQ.setFileExt('.fcm_uniQ')
3908
3909 # Adition of the non-resonant nodes (N+4 matrix)
3910 MatQ.inflateMatrix('both')
3911
3912 # A = -G1 = Im(Mfcm4_22)
3913 A=np.imag(MatQ.M[2, 2])
3914 # B = 2M_12= 2Mfcm4_23
3915 B=2*MatQ.M[2, 3]
3916 # C = - (M_01)^2 / G = (Mfcm4_12)^2 / Im(Mfcm4_11)
3917 C=MatQ.M[1, 2]**2 /np.imag(MatQ.M[1, 1])
3918
3919 # rotation angle alpha from TN.102.3 eq. (10)
3920 alpha=0.5*np.log(1./2/(-2*A-C+2*B)*(-2.*C+4.*np.sqrt(-A**2-C*A+B**2)))
3921
3922 # Aplication of the rotations to make Q uniform
3923 i=2;j=3;
3924 MatQ.rotateMatrix('hyperbolic',i,j,alpha, flagQ=False)
3925 i=5;j=4;
3926 MatQ.rotateMatrix('hyperbolic',i,j,alpha, flagQ=False)
3927
3928 # Node Scaling to make the sum of the imaginary parts of the non-resonant nodes equal to zero, and therefore Q=0 in those nodes
3929 h=-np.imag(MatQ.M[3, 1])/np.imag(MatQ.M[1, 1])
3930 MatQ.scaleNode(1, h, flagQ=False)
3931 h=-np.imag(MatQ.M[4, 6])/np.imag(MatQ.M[6, 6])
3932 MatQ.scaleNode(6, h, flagQ=True)
3933
3934 return MatQ
3935
3936 ##
3937 #
3938 # This function transforms the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix of a filter of order 6 into an N+4 matrix so that all the resonant nodes have equal quality factor Q.
3939 # (Valid for Butterworth, Chevyshev, Quasieliptic and Generalized Chebyshev with symmetric transmission zeros).
3940 # @param MatQfcm2 = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix.
3941 # @return MatQfcm4 = MatrixQ class instance, containing the Folded \f$N+4\f$ Canonical Form (FCM) coupling matrix with uniform Q.
3942 #
3943 # See TN102.3 sec. 3.3.
3944 #
3945 def uniformQFCMOrder6(self, MatQfcm2):
3946 MatQ = MatQfcm2.copy()
3947 MatQ.setName('Folded with uniform Q (N+4)')
3948 MatQ.setFileExt('.fcm_uniQ')
3949
3950 # Obtain the rotation angle which makes Q uniform, with an optimization routine
3951 print "\nUNIFORM Q OPTIMIZATION:"
3952 angRotOpt = scop.fmin(self.uniformQOrder6Optimize, np.array([0,0]), args=(MatQ, 0), xtol=1.e-8, ftol=1.e-8, maxiter=20000, maxfun=20000, full_output=0, disp=True, retall=0, callback=None)
3953 print ""
3954
3955 # Resulting matrix for the optimal angle found
3956 MatQ = self.uniformQOrder6Optimize(angRotOpt, MatQ, 1)
3957
3958 # Conversion to the N+4 matrix taking into account the node scaling done k
3959 MatQ.inflateMatrix('both')
3960
3961 return MatQ
3962
3963
3964 ##
3965 #
3966 # This function performs the hyperbolic rotations in TN102.3 sec.3.3.
3967 # Given two angles angRot it performs the hyperbolic rotations with those angles
3968 # and sees the differences between the quality factors of the different nodes [TN102.3 sec.3.3].
3969 # @param angRot = Vector containing two angles. The first angle corresponds with \f$\alpha\f$ in TN102.3 eq. (14) and the second one with \f$\alpha_{1}\f$
3970 # @param MatQfcm2a = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix.
3971 # @param flag = If False, return goal; if True, return MatQfcm2.
3972 # @return goal = Function we are minimizing. It contains the summation of the absolute values of the diferences in losses in the different nodes.
3973 # @return MatQfcm2 = MatrixQ class instance, containing the transformed Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix.
3974 #
3975 def uniformQOrder6Optimize(self, angRot, MatQfcm2a, flag):
3976 MatQfcm2 = MatQfcm2a.copy()
3977 angle12 = angle65 = angRot[0]
3978 angle23 = angle54 = angRot[1]
3979
3980 # Application of the rotations
3981 i=1; j=2; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle12, flagQ=False)
3982 i=6; j=5; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle65, flagQ=False)
3983 i=2; j=3; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle23, flagQ=False)
3984 i=5; j=4; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle54, flagQ=False)
3985
3986 #Node scaling
3987 k=-np.imag(MatQfcm2.M[2,0])/np.imag(MatQfcm2.M[0,0])
3988 MatQfcm2.scaleNode( 0, k, flagQ=False)
3989 MatQfcm2.scaleNode( 7, k, flagQ=True)
3990
3991 Q = np.array(MatQfcm2.Q)
3992 # Function to minimize
3993 goal = abs(1/Q[0]-1/Q[1]) + abs(1/Q[1]-1/Q[2]) + abs(1/Q[2]-1/Q[3]) + abs(1/Q[3]-1/Q[4]) + abs(1/Q[4]-1/Q[5])
3994
3995 if flag:
3996 return MatQfcm2
3997 else:
3998 return goal
3999
4000
4001 ##
4002 #
4003 # This function transforms the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix of a filter of order 6 into an N+4 matrix so that all the resonant nodes have equal quality factor Q.
4004 # (Valid for Butterworth, Chevyshev, Quasieliptic and Generalized Chebyshev with symmetric transmission zeros).
4005 # @param MatQfcm2 = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix.
4006 # @return MatQfcm4 = MatrixQ class instance, containing the Folded \f$N+4\f$ Canonical Form (FCM) coupling matrix with uniform Q.
4007 #
4008 def uniformQFCMOrder6Case3(self, MatQfcm2):
4009 MatQ = MatQfcm2.copy()
4010 MatQ.setName('Folded with uniform Q (N+4)')
4011 MatQ.setFileExt('.fcm_uniQ')
4012
4013 # Obtain the rotation angle which makes Q uniform, with an optimization routine
4014 print "\nUNIFORM Q OPTIMIZATION:"
4015 angRotOpt = scop.fmin_slsqp(self.uniformQOrder6Case3Optimize, np.array([0,0, 0, 0]) , eqcons=[], f_eqcons=None, ieqcons=[], f_ieqcons=None,bounds = [],
4016 fprime = None, fprime_eqcons=None, fprime_ieqcons=None, args = (MatQ, False), iter = 20000, acc = 1.0E-10, iprint = 1, full_output = 0)
4017 print ""
4018
4019 # Resulting matrix for the optimal angle found
4020 MatQ = self.uniformQOrder6Case3Optimize(angRotOpt, MatQ, True)
4021
4022 # Conversion to the N+4 matrix taking into account the node scaling done k
4023 MatQ.inflateMatrix('both')
4024 return MatQ
4025
4026
4027 ##
4028 #
4029 # This function performs the hyperbolic rotations in TN102.3.
4030 # Given four angles angRot it performs the hyperbolic rotations with those angles.
4031 # and sees the differences between the quality factors of the different nodes [TN102.3].
4032 # @param angRot = Vector containing two angles. The first angle corresponds with \f$\alpha\f$ in TN102.3 eq. (14) and the second one with \f$\alpha_{1}\f$
4033 # @param MatQfcm2a = MatrixQ class instance, containing the Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix.
4034 # @param flag = If False, return goal; if True, return MatQfcm2.
4035 # @return goal = Function we are minimizing. It contains the summation of the absolute values of the diferences in losses in the different nodes.
4036 # @return MatQfcm2 = MatrixQ class instance, containing the transformed Folded \f$N+2\f$ Canonical Form (FCM) coupling matrix.
4037 #
4038 def uniformQOrder6Case3Optimize(self, angRot, MatQfcm2a, flag):
4039 MatQfcm2 = MatQfcm2a.copy()
4040
4041 angle14 = angRot[0];
4042 angle23 = angRot[1];
4043 angle12 = angRot[2];
4044 angle34 = angRot[3];
4045
4046 # Application of the rotations
4047 i=2; j=5; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle14, flagQ=False)
4048 i=3; j=4; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle23, flagQ=False)
4049 i=2; j=3; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle12, flagQ=False)
4050 i=4; j=5; MatQfcm2.rotateMatrix('hyperbolic', i, j, angle34, flagQ=True)
4051
4052 Q = MatQfcm2.Q
4053 # Function to minimize: Q of the four middle resonators, not the 1st and the last ones.
4054 goal = abs(1/Q[1]-1/Q[2]) + abs(1/Q[2]-1/Q[3]) + abs(1/Q[3]-1/Q[4])
4055
4056 if flag:
4057 return MatQfcm2
4058 else:
4059 return goal
4060
4061 # End class: CouplingMatrices
4062 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
4063
4064 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
4065 # Begin: Frequency Transformation classes
4066
4067 ##
4068 #
4069 # Frequency transformation fro band pass to low pass prototype.
4070 #
4071 class FrequencyTransformBP(object):
4072 ##
4073 #
4074 # Create a frequency transformation instance.
4075 #
4076 def __init__(self, P):
4077 self.BW = P.BW
4078 self.f0 = P.f0
4079
4080 ##
4081 # Fractional bandwidth.
4082 self.FBW = P.BW/P.f0
4083
4084
4085 ##
4086 #
4087 # Normalize frequency to lowpass prototype.
4088 # @param freq = Vector with frequency values to convert (in Hz). Numpy array.
4089 # @return Vector with lowpass prototype normalized frequency values. Numpy array.
4090 #
4091 # [Cameron eq. (3.113)]:
4092 #
4093 # \f$ \omega ' = \frac{f_0}{\Delta f} \left( \frac{f}{f_0} - \frac{f_0}{f} \right) \f$
4094 #
4095 def normFreq(self, freq):
4096 # # Obsolete code. It can be removed unless proved otherwise
4097 # if type(freq) != types.ListType: freq = [ freq ]
4098 # ret_freq = [ (self.f0 / self.BW) * (f/self.f0 - self.f0/f) if f is not None else None for f in freq ]
4099 # if len(ret_freq) == 1: ret_freq = ret_freq[0]
4100
4101 return (self.f0 / self.BW) * (freq/self.f0 - self.f0/freq)
4102
4103
4104 ##
4105 #
4106 # Unnormalize frequency from lowpass prototype.
4107 # @param w = Vector with normalized frequency values to convert. Numpy array.
4108 # @return Vector with unnormalized frequency values (in Hz). Numpy array.
4109 #
4110 # \f[
4111 # \omega ' = \frac{f_0}{\Delta f} \left( \frac{f}{f_0} - \frac{f_0}{f} \right)
4112 # \f]
4113 #
4114 # \f[
4115 # f^2 - \omega ' \Delta f f - f_0^2 = 0
4116 # \f]
4117 #
4118 # \f[
4119 # f = \frac {\omega ' \Delta f + \sqrt{( \omega ' \Delta f)^2 + 4 f_0^2 } } {2}
4120 # \f]
4121 # with \f$ \omega'' = \omega' \Delta f / f_0 \f$ :
4122 # \f[
4123 # f = \frac {f_0}{2} \left( \omega'' + \sqrt{ \omega'' \, ^2 + 4 } \right)
4124 # \f]
4125 #
4126 def unormFreq(self, w):
4127 # b = -w*self.BW
4128 # return (-b+np.sqrt(b**2+4*self.f0**2))/2
4129 wp = w*self.FBW
4130 return (self.f0/2) * (wp + np.sqrt(wp**2 +4))
4131
4132
4133 ##
4134 #
4135 # Convert lowpass prototype group delay to unnormalized group delay.
4136 # @param freq = Vector with unnormalized frequency values (in Hz). Numpy array.
4137 # @param groupDelay = Vector with lowpass prototype group delay values to convert. Numpy array.
4138 # @return Vector with unnormalized filter group delay values (in seconds). Numpy array.
4139 #
4140 # [Cameron, eq. (3.90)]:
4141 #
4142 # \f$ \tau_{BPF}= \frac{\omega_0}{\Delta \omega} \left( \frac{1}{\omega_0} - \frac{\omega_0}{\omega^2} \right) \tau ' = \frac{1}{2 \pi \Delta f} \left( 1 + \left( \frac{f_0}{f} \right)^2 \right) \tau ' \f$
4143 #
4144 def unormGroupDelay(self, freq, groupDelay):
4145 return groupDelay * (1 + (self.f0/freq)**2) / (2*pi*self.BW)
4146
4147 ##
4148 #
4149 # Convert unnormalized group delay to lowpass prototype group delay.
4150 # @param groupDelay = Double with unnormalized group delay (in seconds) value to convert.
4151 # @return groupDelay = normalized filter group delay value.
4152 #
4153 # [Cameron, eq. (3.91)]:
4154 #
4155 # \f$ \tau_{0}= \frac{2}{\Delta \omega} \tau_{0} ' \f$
4156 #
4157 def normGroupDelay(self, groupDelay):
4158 return groupDelay * (2*pi*self.BW) / 2
4159
4160 # End: Frequency Transformation classes
4161 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
4162
4163