Microwave filters GUI 2.0.3

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

- File List

libfreefilters.py
Go to the documentation of this file.
1 #/usr/bin/env python
2 # -*- coding: utf-8 -*-
3
4 ## @file
5
6 ##
7 # @package libfreefilters Free Filters Library to compute Butterworth, Chebyshev and Quasieliptic
8 # characteristic polynomials with minimum insertion loss.
9 # @author J.M. Rius, J. Mateu, J.M. Tamayo, C. Collado, A. Padilla and J. O'Callaghan
10 # <p>Dpt. Signal Tehory and Communications, Universitat Politècnica de Catalunya (UPC),
11 # @version 2.0.3
12 # @date June 11, 2013
13 # <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.
14 # Further features were developed under contract UPC-C7767 with Thales Alenia Space España (TAS-E).
15 # Contributions to the definition of the software functionality and testing have been made by Christoph Ernst, Mónica Martínez Mendoza and other
16 # ESA-ESTEC personnel, and Santiago Sobrino and Luis Roglá from TAS-E.
17 # <p>Copyright: &copy; Universitat Politècnica de Catalunya (UPC) 2009. <br>Contact: lossyfilters@tsc.upc.edu
18 # <p>This program is free software; you can redistribute it and/or modify it under the terms of
19 # the GNU General Public License as published by the Free Software Foundation; either
20 # version 3 of the License, or (at your option) any later version.
21 # <p>This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
22 # without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
23 # See the GNU General Public License for more details.
24 # <p>You should have received a copy of the GNU General Public License along with this program
25 # in file "LICENSE.GPL3"; if not, download it from http://www.gnu.org/licenses/gpl-3.0.html .
26 # <p>Additional permission under GNU GPL version 3 section 7:
27 # <p>If you modify this Program, or any covered work, by linking or combining it with the "Extra Filters Library" (libextrafilters),
28 # the "Lossy Filters Library" (liblossyfilters) or the "License Check Library" (libchecklicense), or modified versions of that libraries,
29 # containing parts covered by the terms described in files LICENSE.LIBEXTRAFILTERS and LICENSE.LIBLOSSYFILTERS,
30 # the licensors of this Program grant you additional permission to convey the resulting work.
31 #
32
33 from math import log10, pi, sin, cos, sqrt, exp, sinh, cosh
34 import types
35 import scipy as sc
36 import numpy as np
37 import os.path
38 import platform
39 from PyQt4.QtGui import QMessageBox
40
41 from libcommonfunc import pkgNameVersion, complexStr, myPrint, printHeader, synthError, askQuestion, readVerbose, warningMsg
42
43 printHeader('Successful import:' , __doc__, importModule = True)
44
45 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
46 # Import libextrafilters AT THE END of the file, to avoid CIRCULAR IMPORT #
47 # Import mwfiltersgui AT THE END of the file, to avoid CIRCULAR IMPORT #
48 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
49
50 ##
51 # Module global variable equal either to None or to the GUI QMainWindow class instance, containing the application main window.
52 # It is equal to None (default) when the calling program is the console interface (mwfiltersgui.py --nogui).
53 # When creating the main window, the GUI executes the libfreefilters.setGlobalVariablesLibFreeFilters() function in order to set this variable.
54 # The GUI methods do not use this variable. They access the mainWindow class through ParamEditDlg.mainWindow and MainWindow 'self' attributes.
55 mainWindow = None
56
57 ##
58 #
59 # Function to allow the GUI set values to the global variables of libcomonfunc module.
60 # This function is used in order to avoid circular imports.
61 # @param variablesDict: Variable keyword argument list. Python creates a dictionary.
62 #
63 def setGlobalVariablesLibFreeFilters(**variablesDict):
64 for key, value in variablesDict.items():
65 globals()[key] = value
66
67 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
68 # Begin class: CharPolys #
69
70 ##
71 #
72 # Characteristic polynomials class.
73 #
74 # Characteristic polynomials and constants are defined as [TN 101 eq. (1)] :
75 # \f[
76 # S_{11}(s) = \left[ \frac{F(s)}{\epsilon_r E(s)} \right],
77 # S_{12}(s) = S_{21} = \left[ \frac{P(s)}{\epsilon E(s)} \right],
78 # S_{22}(s) = \left[ \frac{(-1)^N F(s)^*}{\epsilon_r E(s)} \right]
79 # \f]
80 # where \f$ F(s)^* = F^*(-s) \f$ is the 'Paraconjugate' polynomial of \f$ F(s) \f$ [Cameron pp. 208 footnote]
81 # and is computed using polyConj function.
82 #
83 class CharPolys(object):
84 ##
85 #
86 # Compute characteristic polynomials and characteristic constants.
87 # @param P = Parameter class instance containing filter, synthesis and sweep parameters.
88 # @param FT = Frequency transformation class instance.
89 # @param symmetrizeZeros = Automatically symmetrize Generalized Chebyshev zeros in Hz in order to compute a folded matrix with uniform Q in resonators.
90 # Equivalent to answering "Yes" in the GUI when it asks the user if he wants to symmetrize Generalized Chebyshev zeros in Hz.
91 # @return CharPolys class instance.
92 #
93 def __init__(self, P, FT, symmetrizeZeros = None):
94
95 ##
96 # Characteristic polynomial \f$ P(s) \f$
97 self.Ps = None
98
99 ##
100 # Characteristic polynomial \f$ E(s) \f$
101 self.Es = None
102
103 ##
104 # Characteristic polynomial \f$ F(s) \f$
105 self.Fs = None
106
107 ##
108 # Characteristic polynomial \f$ F_{22}(s) \f$
109 self.F22s = None
110
111 ##
112 # Characteristic constant \f$ \epsilon \f$ in the denominator of \f$ S_{12} \f$ and \f$ S_{21} \f$
113 self.eps = None
114
115 ##
116 # Characteristic constant \f$ \epsilon_R \f$ in the denominator of \f$ S_{11} \f$
117 self.epsR = None
118
119 ##
120 # Roundoff error estimation in the computation of lossless characteristic polynomials,
121 # given by the rootsErr return value in CharPolys.factorization method.
122 self.rootsErr = None
123
124 ##
125 # Value 'a' of the (s-a)/(s+a) zero/pole added in lossy filters synthesis case 3.
126 self.nuqDelta = None
127
128 # Filter order
129 N = P.N
130
131 if N < 2: raise synthError, "Filter order must be equal or greater than 2"
132
133 # Set flags for computation of Uniform Q in special topologies of TN 102.3 (cases 1, 2 or 3 )
134 # Uniform Q optimization in cases 1, 2 or 3 can be performed only if TZs are symmetric in normalized frequency
135 self.flagExtraFilter = self.flagSymTZ = self.flagSymTZnorm = self.flagSymTZunorm = self.flagComplexTZ = self.flagCase1 = self.flagCase2 = self.flagCase3 = False
136
137 # Standard filters ( u'Butterworth', u'Chebyshev', u'Quasieliptic' ) have no TZs or symmetric TZs in normalized frequency
138 if P.filterTransferFunc.title() in ( u'Butterworth', u'Chebyshev', u'Quasieliptic' ):
139 self.flagSymTZ = True
140 else: # Only for filters in the 'Extra Filter Library', we have to find out if TZs are symmetric or not
141 self.flagExtraFilter = True
142
143 # Be sure that lists of zeros are really lists
144 if type(P.transImagZeros) != types.ListType: P.transImagZeros = [ P.transImagZeros ]
145 if type(P.transComplexZeros) != types.ListType: P.transComplexZeros = [ P.transComplexZeros ]
146 if type(P.LPcomplexZeros) != types.ListType: P.LPcomplexZeros = [ P.LPcomplexZeros ]
147 if type(P.LPrealZeros) != types.ListType: P.LPrealZeros = [ P.LPrealZeros ]
148
149 # Check that all TZs are outside the passband
150 if np.any( np.abs(FT.normFreq(np.array( P.transImagZeros ))) < 1):
151 raise synthError, 'There are transmission zeros in the passband.'
152
153 # Check if TZs are symmetric in normalized and in unnormalized frequency axis.
154 self.flagSymTZunorm = self.checkUnormSymTZ(P.transImagZeros, FT)
155 self.flagSymTZnorm = self.checkNormSymTZ(FT.normFreq(np.array(P.transImagZeros)))
156 self.flagSymTZ = self.flagSymTZnorm or self.flagSymTZunorm
157
158 # Set complex transmission zeros flag, must be False for 'Butterworth', 'Chebyshev', 'Quasieliptic'
159 if P.genChebyLP or len(P.transComplexZeros) > 0: self.flagComplexTZ = True
160
161
162 # Set case1, case2 or case3 flag if transmission zeros are symmetric and not complex
163 if self.flagSymTZ and not self.flagComplexTZ: # No transmission zeros or symmetric transmission zeros
164 if (N == 4) and (P.synthTech.title() == u'Prescribed Insertion Loss') and (P.nuqTech == u'kS21+kS11') and (P.nuqK11c1 == P.nuqK21c1 and P.nuqK11c1 == P.nuqK22c1):
165 self.flagCase1 = True
166 elif (N == 6) and (P.synthTech.title() == u'Prescribed Insertion Loss') and (P.nuqTech == u'kS21+kS11') and (P.nuqK11c1 == P.nuqK21c1 and P.nuqK11c1 == P.nuqK22c1):
167 self.flagCase2 = True
168 elif (N == 6) and (P.synthTech.title() == u'Prescribed Insertion Loss') and (P.nuqTech == u'kS21+pole/zero'):
169 self.flagCase3 = True
170
171 # If TZs are symmetric, but not in normalized frequency, we have to symmetrize them, but first ask the user if he wants to
172 if self.flagExtraFilter and any([self.flagCase1, self.flagCase2, self.flagCase3]) and not self.flagSymTZnorm:
173 if symmetrizeZeros is None: # The execution comes from the GUI and symmetrizeZeros command line flag is not set.
174 reply = askQuestion(
175 """<center><b>Asymmetric transmission zeros:</b></center><p>
176 In order to compute a folded coupling matrix with uniform Q of resonators, the software needs that transmission zeros are symmetric in the normalized frequency axis.<p>
177 Transmission zeros are symmetric with respect to center frequency in Hz, but not in normalized frequency.<p>
178 <center><font color="red"><b>Do you want to move transmission zeros?</b></font></center>""",
179 QMessageBox.Yes | QMessageBox.No )
180
181 symmetrizeZeros = (reply == QMessageBox.Yes)
182
183 if symmetrizeZeros:
184 P.transImagZeros = self.symmetrizeTZ(P.transImagZeros, FT)
185 else: # The user does not want to symmetrize zeros, so we have to disable Uniform Q optimization
186 self.flagCase1 = self.flagCase2 = self.flagCase3 = False
187
188 # Set the roundoff error in the computation with the current precision
189 self.rootsErr = sc.arcsin(sin(pi))
190
191 if P.filterTransferFunc.title() == u'Butterworth': Ps, Es, Fs, eps, epsR, rootsErr = butterworthFilter(P, FT, symmetrizeZeros)
192 elif P.filterTransferFunc.title() == u'Chebyshev': Ps, Es, Fs, eps, epsR, rootsErr = chebyshevFilter(P, FT, symmetrizeZeros)
193 elif P.filterTransferFunc.title() == u'Quasieliptic': Ps, Es, Fs, eps, epsR, rootsErr = quasiElipticFilter(P, FT, symmetrizeZeros)
194
195 elif P.filterTransferFunc.title() == u'Generalized Chebyshev':
196 if importedExtraFilters is False:
197 if existingExtraFilters is False:
198 raise synthError, """Generalized Chebyshev filter is in the Library of Extra Filters,<br>
199 which is not installed in this system.
200 <p>Contact lossyfilters@tsc.upc.edu to get the Library of Extra Filters.
201 """
202 else:
203 raise synthError, """Importing the 'Library of Extra Filters' failed. <br>
204 You have a libextrafilters.pyc file, but it cannot be imported. <br>
205 Possibly libextrafilters.pyc was compiled using a different Python version than yours.
206 <p>Your installed Python version is: %s on %s
207 <p>Read the README.html documentation or contact lossyfilters@tsc.upc.edu for assistence.
208 """ % (platform.python_version(), platform.system())
209 elif importedCheckLicense is False:
210 if existingCheckLicense is False:
211 raise synthError, """There is no check_license.pyc file in your system.
212 It should have been installed together with libextrafilters.pyc file,
213 which actually is in your system.
214 Please follow the 'Extra Filters Library' installation instructions.
215 """
216 else:
217 raise synthError, """Importing the 'License Check Library' failed. <br>
218 You have a libchecklicense.pyc file, but it cannot be imported. <br>
219 Possibly libchecklicense.pyc was compiled using a different Python version than yours.
220 <p>Your installed Python version is: %s on %s
221 <p>Read the README.html documentation or contact lossyfilters@tsc.upc.edu for assistence.
222 """ % (platform.python_version(), platform.system())
223 else:
224 Ps, Es, Fs, eps, epsR, rootsErr = generalizedChebyshevFilter(P, FT, symmetrizeZeros)
225
226 else: raise synthError, u"Filter topology type '%s' not available." % P.filterTransferFunc
227
228 # Update roundoff error
229 self.rootsErr = max(rootsErr, self.rootsErr)
230
231 # Compute F22(s) according to [Cameron eq. (6.15)]
232 F22s = (-1)**Fs.o * polyConj(Fs)
233
234 # Copy local variables to class attributes
235 self.Ps, self.Es, self.Fs, self.F22s, self.eps, self.epsR = Ps, Es, Fs, F22s, eps, epsR
236
237 ##
238 #
239 # A list of transmission zeros in unnormalized frequency (Hz) is symmetrized in the normalized frequency axis.
240 # The positive zeros are averaged with the corresponding negative zeros in the normalized frequency axis.
241 # @param TZ = List of transmission zeros. Unnormalized frequency (Hz).
242 # @param FT = Frequency transform class instance.
243 # @return TZsym = List of transmission zeros, symmetrized n the normalized frequency axis. Unnormalized frequency (Hz).
244 #
245 def symmetrizeTZ(self, TZ, FT):
246 tol = 0
247 TZ = np.array(TZ)
248 TZnorm = np.array(FT.normFreq(TZ))
249 TZpos = np.sort(TZnorm[np.nonzero(TZnorm>tol)])
250 TZneg = np.sort(-TZnorm[np.nonzero(TZnorm<-tol)])
251 TZnull = len(np.nonzero(np.abs(TZnorm)<=tol)[0])
252 TZsymNorm = (TZpos+TZneg)/2.
253 TZsymNorm = np.concatenate((np.zeros(TZnull), TZsymNorm, -TZsymNorm))
254 TZsym = [ f for f in FT.unormFreq(TZsymNorm) ]
255 return TZsym
256
257 ##
258 #
259 # Check if list of transmission zeros (in Hz) is symmetric with respect to center frequency f0.
260 # @param TZ = List of transmission zeros in Hz.
261 # @param FT = Frequency transform class instance. It will be used to get f0.
262 # @return symmetricFlag = True is zeros are symmetric with relative error < 1e-3, False otherwise.
263 #
264 def checkUnormSymTZ(self, TZ, FT):
265 if len(TZ) == 0: return True
266 tol = 0
267 TZ = np.array(TZ)
268 TZpos = np.sort(TZ[np.nonzero(TZ - FT.f0 > tol)])
269 TZneg = np.sort(TZ[np.nonzero(TZ - FT.f0 < -tol)])
270 if len(TZpos)==len(TZneg):
271 if np.max(np.abs((TZpos+TZneg)/(2*FT.f0)-1)) < 1e-3:
272 return True
273 return False
274
275
276 ##
277 #
278 # Check if list of normalized transmission zeros is symmetric with respect to zero.
279 # @param TZ = List of normalized transmission zeros.
280 # @return symmetricFlag = True is zeros are symmetric with relative error < 1e-3, False otherwise.
281 #
282 def checkNormSymTZ(self, TZ):
283 if len(TZ) == 0: return True
284 tol = 0
285 TZ = np.array(TZ)
286 TZpos =np.sort(TZ[np.nonzero(TZ > tol)])
287 TZneg = np.sort(TZ[np.nonzero(TZ < -tol)])
288 if len(TZpos)==len(TZneg) and np.max(np.abs(TZpos+TZneg)) < 1e-3:
289 return True
290 else:
291 return False
292
293
294 ##
295 #
296 # Return string with polynomial coefficients separated by commas.
297 # Rounds polynomial coefficients to given precision using complexStr() function.
298 # @param pol = Polynomial (numpy.poly1d object).
299 # @param prec = Number of significant digits (int).
300 # @return stcoef = String representation of pol.
301 #
302 def polyStr(self, pol, prec):
303 stcoef = ''
304 for k in range(pol.order+1):
305 #re, im = real(pol.coef[k]), imag(pol.coef[k])
306 #stre = '%.*g' % (prec, re) if abs(re) > 10**-prec else '0'
307 #stim = '% +.*gj' % (prec, im) if abs(im) > 10**-prec else ''
308 #stcoef += stre + stim + ', '
309 stcoef += complexStr(pol.coef[k], prec) + ', '
310 return stcoef[0:-2]
311
312
313 ##
314 #
315 # Save Characteristic Polynomials and constants in output file.
316 # Rounds polynomial coefficients to given precision using complexStr() function.
317 # Polynomials are of type numpy.poly1d.
318 # @param P = Parameter class instance.
319 # @param prec = Number of correct significant digits to print.
320 #
321 def saveCharPolys(self, P, prec):
322 fileName = P.outDirName + '.pol'
323 fout = open(fileName, 'w')
324 fout.write('# %s; version %s\n' % pkgNameVersion(__doc__) )
325 fout.write('# Characteristic polynomials file.\n')
326 fout.write('# Polynomial coefficients are ordered from highest to lowest degree.\n\n')
327 fout.write( 'F(s) = %s\n' % self.polyStr(self.Fs, prec))
328 fout.write( 'P(s) = %s\n' % self.polyStr(self.Ps, prec))
329 fout.write( 'E(s) = %s\n' % self.polyStr(self.Es, prec))
330
331 fout.write('\n# Constant in the denominator of S12 and S21\n')
332 fout.write('eps = %s\n' % complexStr(self.eps, prec) )
333
334 fout.write('\n# Constant in the denominator of S11\n')
335 fout.write('epsR = %s\n' % complexStr(self.epsR, prec) )
336 fout.close()
337 myPrint("Characteristic Polynomials saved to file: %s" % fileName )
338
339 # End class: CharPolys
340 #!!!!!!!!!!!!!!!!!!!!!!!!!!!
341
342
343 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
344 # Begin auxiliary functions
345
346 ##
347 #
348 # Returns \f$ P_s(s) \f$ for the input polynomial \f$ P_w(\omega) \f$, with \f$ s=j\omega \f$ or \f$ \omega=-js \f$.
349 # If \f$ N \f$ is the polynomial order, \f$P_w(\omega)\f$ and \f$ P_s(s) \f$ can be expressed as:
350 # \f[
351 # P_w(w) = \sum_{n=0}^{N}c_{n}\omega^{n}
352 # \f]
353 # \f[
354 # P_s(s) = [P_w(\omega)]|_{\omega=-js} = \sum_{n=0}^{N}(-j)^{n}c_{n}s^{n}
355 # \f]
356 # @param Pw = Input polynomial \f$ P_w(\omega) \f$.
357 # @return Ps = \f$ P_s(s) \f$ .
358 #
359 def polyOmega(Pw):
360 N = Pw.o # Order of the polynomial
361 Ps = np.poly1d(Pw.coef * pow(-1j, np.arange(N+1.)[::-1])) # Computation of the conjugate
362 return Ps
363
364
365 ##
366 #
367 # Compute Chebyshev polynomials from order 0 to N.
368 # @param N = Maximum order of polynomials (int)
369 # @return Tn = List containing Chebyshev polynomials. List elements of type numpy.poly1d.
370 #
371 # \section alfPolyCheby Algorithm
372 #
373 # The recursion formula for Chebyshev polynomials
374 # \f[
375 # T_{n}(\omega) = 2 \omega T_{n-1}(\omega) - T_{n-2}(\omega)
376 # \f]
377 # is used to find \f$ T_N \f$ after setting \f$ T_0 = 1 \f$ and \f$ T_1 = \omega \f$.
378 #
379 def polyCheby(N):
380 Tn = [ np.poly1d([1]), np.poly1d([1, 0]) ]
381 for k in range(2, N+1): Tn.append( 2*Tn[1]*Tn[k-1] - Tn[k-2] )
382 return Tn
383
384
385 ##
386 #
387 # Returns the 'Paraconjugate' \f$ P_c(s) \f$ of the input polynomial \f$ P(s) \f$, with complex variable restricted to the imaginary axis, \f$ s=j\omega \f$, and therefore \f$ s^{*} = -s \f$.
388 # If \f$ N \f$ is the polynomial order, \f$P\f$ can be expressed as:
389 # \f[
390 # P(s) = \sum_{n=0}^{N} c_{n}s^{n}
391 # \f]
392 # We define the paraconjugate \f$P_c(s)\f$ as [Cameron pp. 208, footnote].:
393 # \f[
394 # P_c(s) = P(s)^* = P^*(s^*) = P^{*}(-s) = \sum_{n=0}^{N}(-1)^{n}c_{n}^*s^{n}
395 # \f]
396 # When the frequency is evaluated for complex \f$ s \f$ not in the imaginary axis, \f$ s = \sigma + j\omega \f$, you cannot compute \f$ S_{22} \f$ using
397 # the 'paraconjugate' polynomial of \f$ F(s) \f$, and must complex-conjugate the array of values of \f$ F(s) \f$ for each \f$ s \f$ .
398 # @param P = Input polynomial \f$ P(s) \f$.
399 # @return Pc = Paraconjugate polynomial \f$P_c(s)\f$.
400 #
401 def polyConj(P):
402 N = P.o # Order of the polynomial
403 Pc = np.poly1d(np.conj(P)*pow(-1, np.arange(N+1.)[::-1])) # Computation of the conjugate
404 return Pc
405
406
407 ##
408 #
409 # This function finds the maximum of the function \f$ g(s)=|\frac{NUM(s)}{\epsilon DEN(s)}| \f$ in the imaginary axis. Usually \f$ NUM = Ps,Fs \f$ and \f$ DEN=Es \f$.
410 # Polynomials are of type numpy.poly1d.
411 # @param NUM = Numerator polynomial of the function to maximize.
412 # @param DEN = Denominator polynomial of the function to maximize.
413 # @return eps = Maximum value of the function to maximize. It accomplishes that the maximum value of \f$ g(s)=|\frac{NUM(s)}{\epsilon DEN(s)}| \f$ in the imaginary axis is equal to 0 dB.
414 # To do the maximization we procede as follows:
415 # \f$ SQNUM(s) =|NUM(s)|^{2} = NUM(s)NUM^{*}(-s) \f$ for \f$ s=j\omega \f$.
416 # \f$ SQDEN(s)=|DEN(s)|^{2} = DEN(s)DEN^{*}(-s) \f$ for \f$ s=j\omega \f$.
417 # The function f to maximize is:
418 # \f[ f(s)=\frac{SQNUM(s)}{SQDEN(s)}=\frac{|NUM(s)|^{2}}{|DEN(s)|^{2}} \f]
419 # \f[ f'(s)=\frac{SQNUM'(s)SQDEN(s)-SQNUM(s)SQDEN'(s)}{SQDEN(s)^{2}}=0 \f]
420 # The maximum of f is obtained in one of the roots of the numerator of \f$ f'(s) \f$.
421 #
422 def findEps(NUM, DEN):
423 # \f$ SQNUM(s) =|NUM(s)|^{2} = NUM(s)NUM^{*}(-s) \f$ for \f$ s=j\omega \f$.
424 SQNUM = NUM * polyConj(NUM)
425 # \f$ SQDEN(s)=|DEN(s)|^{2} = DEN(s)DEN^{*}(-s) \f$ for \f$ s=j\omega \f$.
426 SQDEN = DEN * polyConj(DEN)
427 # The function f to maximize is:
428 # \f[ f(s)=\frac{SQNUM(s)}{SQDEN(s)}=\frac{|NUM(s)|^{2}}{|DEN(s)|^{2}} \f]
429 # \f[ f'(s)=\frac{SQNUM'(s)SQDEN(s)-SQNUM(s)SQDEN'(s)}{SQDEN(s)^{2}}=0 \f]
430 # The maximum of f is obtained in one of the roots of the numerator of \f$ f'(s) \f$.
431 DERnum = SQNUM.deriv() * SQDEN - SQNUM * SQDEN.deriv()
432 jOm = DERnum.r
433 POL=abs(NUM(1j*np.imag(jOm))/DEN(1j*np.imag(jOm)))
434 eps = max(POL)
435 return eps
436
437
438 ##
439 #
440 # Report errors in the roots position for a polynomial of the form SQF(s) = F(s)F(-s) with real coefficients.
441 # The roots of the polynomial SQF(s) must be symmetric about imaginary axis. Note that for pure imaginary roots, they must be of even multiplicity.
442 #
443 # The error is computed as the maximum distance between the roots theoretical and actual positions.
444 # @param roots_SQF = Array of roots. Type numpy.array.
445 # @return maxErr = Maximum distance between the roots theoretical and actual positions.
446 #
447 def rootErrors(roots_SQF):
448 N = len(roots_SQF)
449 maxErr = 0
450 for k in range(N):
451 #errConj = min( abs(roots_SQF - np.conj(roots_SQF[k])) ) # Error in the position of the conjugate
452 #errOp = min( abs(roots_SQF + roots_SQF[k]) ) # Error in the position of the opposite
453 #errOpConj = min( abs(roots_SQF + np.conj(roots_SQF[k])) ) # Error in the position of the conjugate opposite
454 #maxErr = max([ maxErr, errConj, errOp, errOpConj ])
455 errSym = min( np.real ( roots_SQF + roots_SQF[k] ))
456 maxErr = max( maxErr, errSym )
457 return maxErr
458
459
460 ##
461 #
462 # Factorization of polynomial SQF(s) of order 2N as a product of two polynomials, each one of order N.
463 # Polynomials are of type numpy.poly1d.
464 # @param SQF = Polynomial to factorize as a product of two.
465 # @param method = Either an integer value between 1 and 4 describing the method of choosing which roots of the polynomial SQF go in each factor, or a list of Booleans.
466 # <br>See documentation below.
467 # @return Fs, F22s = The two polynomials in which the input polynomial SQF is factorized.
468 # Fs and F22s are normalized to highest order coefficient equal to 1.
469 # F22s is formed from the roots discarded in Fs, and should be equal to
470 # \f$ (-1)^N F(s)^* = (-1)^N F^*(-s) \f$ [Cameron pp. 208,212],
471 # where \f$ F(s)^* \f$ is the 'Paraconjugate' polynomial of \f$ F(s) \f$ (see function CharPolys.polyConj).
472 # @return rootsErr = Maximum distance between the theoretical position of SQF roots and the actual position.
473 # Real or imaginary parts with abs() smaller than tolerance will be considered as a zero.
474 # A pair of roots separated less than this tolerance in s-plane will be cosidered as a double root.
475 #
476 # \section fact Factorization
477 #
478 # The input polynomial is factorized as:
479 # \f[
480 # SQF(s) = F(s)F_{22}(s)
481 # \f]
482 # so that when s is restricted to the imaginary axis, \f$ s=j\omega \f$, it accomplishes the following:
483 # \f[
484 # |SQF(j\omega)| = F(j\omega)F^{*}(j\omega),
485 # |SQF(j\omega)| = F_{22}(j\omega)F_{22}^{*}(j\omega), \forall \omega \in \Re.
486 # \f]
487 # The roots of the polynomial SQF(s) must be symmetric about imaginary axis. Note that for pure imaginary roots, they must be of even multiplicity.
488 # This restriction is equivalent to the fact that SQF(s) can be factorized as:
489 # \f[
490 # SQF(s) = F(s)F(-s)
491 # \f]
492 #
493 # There are different methods to choose the roots of F(s) from the roots of SQF(s): <br>
494 # <ol>
495 # <li>If 'method' argument is an integer:
496 # <ul>
497 # <li>method == 1 -> All zeros in the left s-plane (Hurwitz).
498 # <li>method == 2 -> All zeros in the right s-plane.
499 # <li>method == 3 -> The roots in right and left s-planes are ordered by decreasing imaginary part and alternatively we choose for F(s) a root in the right or in the left s-plane.
500 # <li>method == 4 -> The complementary of 3.
501 # </ul>
502 # <li>If 'method' argument is a list of Booleans:
503 # <ul>
504 # <li>If method[n] is True, zeros_in_left_s-plane[n] go to F(s) and the symmetrical zeros_in_right_s-plane[n] go to F(-s).
505 # <li>If method[n] is False, zeros_in_left_s-plane[n] go to F(-s) and the symmetrical zeros_in_right_s-plane[n] go to F(s).
506 # <li>The zeros zero_in_right_s-plane and zero_in_left_s-plane are ordered by decreasing imaginary part.
507 # </ul>
508 # </ol>
509 # F22(s) is composed of the remaining roots from SQF(s) to form the complementary function to F(s). Method 1 is complementary to method 2 and 3 to 4:
510 # \f[
511 # F^{1}(s) = F^{2}_{22}(s), F^{2}(s) = F^{1}_{22}(s),
512 # F^{3}(s) = F^{4}_{22}(s), F^{4}(s) = F^{3}_{22}(s).
513 # \f]
514 #
515 def factorization(SQF, method):
516
517 N = SQF.o # Degree of the polynomial SQF to factorize.
518 roots_SQF = SQF.roots # Roots of the polynomial SQF to factorize.
519
520 # Assuming quad symmetry of complex roots in SQF, test roundoff error in the positiom of roots
521 rootsErr = rootErrors(roots_SQF)
522
523 # Separation into different kind of roots.
524 roots_PRP = roots_SQF[np.nonzero( np.real(roots_SQF) > rootsErr) ] # Complex roots (not purely imaginary) with positive real part.
525 roots_NRP = roots_SQF[np.nonzero( np.real(roots_SQF) < -rootsErr) ] # Complex roots (not purely imaginary) with negative real part.
526 roots_ZRP = roots_SQF[np.nonzero( np.abs(np.real(roots_SQF)) <= rootsErr) ] # Purely imaginary roots with zero real part (they are double roots)
527 roots_ZPRP = vecSort(roots_ZRP)[::2].copy() # Purely imaginary roots with zero-positive real part (simple roots)
528 roots_ZNRP = vecSort(roots_ZRP)[1::2].copy() # Purely imaginary roots with zero-negative real part (simple roots)
529
530 # Remove round off errors
531 roots_ZPRP = 1j*np.imag(roots_ZPRP)
532 roots_ZNRP = 1j*np.imag(roots_ZNRP)
533
534 # Join vectors of zeros
535 roots_PRP = vecSort(np.concatenate( ( roots_PRP, roots_ZPRP ), 1))
536 roots_NRP = vecSort(np.concatenate( ( roots_NRP, roots_ZNRP ), 1))
537
538 if readVerbose():
539 myPrint('\nFactorization of F(s)F(-s)):')
540 myPrint('Zeros in left-s plane =' + str( [ complexStr(r, 3) for r in roots_NRP ] ) )
541 myPrint('Zeros in right-s plane =' + str( [ complexStr(r, 3) for r in roots_PRP ] ) )
542
543 # Open GUI for selecting zeros of F(s), only if method is already a list of bools
544 if mainWindow is not None and type(method).__name__ == 'list' and type(method[0]).__name__ == 'bool':
545 PZdlg = PredisZeorsDlg(roots_NRP, method, mainWindow)
546 if PZdlg.exec_():
547 pass
548 else:
549 pass
550 del PZdlg
551
552 # Compute the roots for F and F22 depending on the chosen method.
553 if type(method).__name__ == 'int':
554 if method == 1: # F with all roots with negative real part.
555 roots_F = vecSort(roots_NRP.copy())
556 roots_F22 = vecSort(roots_PRP.copy())
557 elif method == 2: # F with all roots with positive real part.
558 roots_F = vecSort(roots_PRP.copy())
559 roots_F22 = vecSort(roots_NRP.copy())
560 elif method == 3:
561 roots_F = vecSort(np.concatenate( (roots_PRP[::2], roots_NRP[1::2]), 1))
562 roots_F22 = vecSort(np.concatenate( (roots_PRP[1::2], roots_NRP[::2]), 1))
563 elif method == 4:
564 roots_F = vecSort(np.concatenate( (roots_PRP[1::2], roots_NRP[::2]), 1))
565 roots_F22 = vecSort(np.concatenate( (roots_PRP[::2], roots_NRP[1::2]), 1))
566 else:
567 raise synthError, u"Factorization of F(s)F(-s): Factorization method is not an integer in the range [1,4]."
568
569 elif type(method).__name__ == 'list' and type(method[0]).__name__ == 'bool':
570 if len(method) != len(roots_PRP): raise synthError, u"Length of Factorization method boolean list is different than filter order."
571 aMethod = np.asarray(method) # Convert method to numpy array, for indexing
572 roots_F = vecSort(np.concatenate( (roots_PRP[ ~aMethod], roots_NRP[ aMethod]), 1))
573 roots_F22 = vecSort(np.concatenate( (roots_PRP[ aMethod], roots_NRP[~aMethod]), 1))
574
575 else: raise synthError, u"Factorization of F(s)F(-s): Factorization method is not an integer in the range [1,4] or a list of Booleans."
576
577 if readVerbose():
578 myPrint('Factorization method = ' + str(method))
579 myPrint('F(s) =' + str( [ complexStr(r, 3) for r in roots_F ] ) )
580 myPrint('F(-s) =' + str( [ complexStr(r, 3) for r in roots_F22 ] ) )
581 myPrint('')
582
583 # Compute polynomials from their roots
584 Fs = np.poly1d(roots_F, r=1)
585 F22s = np.poly1d(roots_F22, r=1)
586
587 # Normalize polynomials to have highest order coef equal to 1
588 Fs = Fs / Fs.coef[0]
589 F22s =F22s / F22s.coef[0]
590
591 # Return the two polynomials taking into account that the first coeffcicient of SQF could be different from one
592 return Fs, F22s , rootsErr
593
594
595 ##
596 #
597 # Sort a vector of complex numbers so that their imaginary parts is in decreasing order.
598 # @param vec = complex vector to sort.
599 # @return vecsort = sorted vector.
600 #
601 def vecSort(vec):
602 # Index order for the imaginary parts in decreasing order
603 ind = np.argsort(np.imag(vec))[: : -1]
604 # Reorder the input vector vec
605 vecsort = vec[ind]
606 return vecsort
607
608
609 # End: auxiliary functions
610 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
611
612
613 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!
614 # Begin filter functions
615
616 ##
617 #
618 # Compute characteristic polynomials and characteristic constants of a Butterwoth filter.
619 # @param P = Parameter class instance containing filter, synthesis and sweep parameters.
620 # @param FT = Frequency transformation class instance.
621 # @param symmetrizeZeros = This parameter is ignored in this class of filter.
622 # @return Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$ over characteristic constant \f$\epsilon\f$ (\f$ P(s)/\epsilon \f$).
623 # @return Es = Characteristic polynomial of the denominator of \f$S_{21}\f$ and \f$S_{11}\f$.
624 # @return Fs = Caracteristic polynomial of the numerator of \f$S_{11}\f$ considering a lossy system.
625 # @return eps = Constant in the denominator of \f$S_{21}\f$.
626 # @return epsR = Constant in the denominator of \f$S_{11}\f$.
627 # @return rootsErr = Maximum distance between the theoretical position of SQF roots and the actual position. Equal to zero for this kind of filter.
628 #
629 # \section butterCharPolys Algorithm
630 # [TN 102.1 sec. 7.2].
631 # <ol>
632 # <li> The roots of \f$ E(s)E(-s) \f$ can be found as [TN 102.1 eq. (28)]:
633 # \f[
634 # s_k = \left\{ \begin{array}{ll}
635 # \exp \left( \frac{j\pi(2k-1)}{2N} \right) & \text{$N$ odd} \\
636 # \exp \left( \frac{j\pi 2k }{2N} \right) & \text{$N$ even}
637 # \end{array} \right.
638 # \f]
639 # where \f$ k = 1 \ldots 2N \f$ and \f$ N \f$ is the filter order.
640 #
641 # The coeffcicients of \f$ E(s) \f$ are obtained from the roots of \f$ E(s)E(-s) \f$ that lie in the LHS of the s-plane using the numpy.poly1d methods.
642 #
643 # <li> \f$ P(s) = 1 \f$
644 #
645 # <li> \f$ F(s) = s^N \f$
646 #
647 # <li> \f$ \epsilon = 1 \f$
648 #
649 # <li> \f$ \epsilon_R = 1 \f$
650 # </ol>
651 #
652 def butterworthFilter(P, FT, symmetrizeZeros = None):
653
654 N = P.N
655
656 if N>=10: warningMsg(u'<b>WARNING:</b> <p>Butterworth filter coupling matrices are not accurate for order 10 or more.')
657
658 # Roots of E(S) [TN 102.1 eq. (28)]
659 k = np.arange(1, 2*N+1)
660 if N%2 == 0: # N even
661 rootsEs = np.exp(1j * pi * (2*k-1) / (2*N) )
662 else: # N odd
663 rootsEs = np.exp(1j * pi * k / N)
664 rootsEs = rootsEs[np.where(np.real(rootsEs) < 0)] # Pick up roots in the left plane
665
666 #Es = np.poly1d( poly(rootsEs), variable='s' )
667 Es = np.poly1d( rootsEs, r=1, variable='s' )
668 Es = Es/Es.coef[0] # Normalize Es
669
670 # P(s): TN 102.1 sec. 7.2
671 Ps = np.poly1d(1., variable='s' )
672 # Orthogonailty law: [David's PFC pp. 11-12], [Cameron pp. 213]
673 if (N - Ps.o)% 2 == 0: Ps *= 1j
674
675 # F(s): TN 102.1 sec. 7.2
676 Fs = np.poly1d(np.concatenate(( np.array([1]), np.zeros(N) )), variable='s' )
677 eps = 1.
678 epsR = 1.
679
680 rootsErr = 0
681
682 return Ps, Es, Fs, eps, epsR , rootsErr
683
684
685 ##
686 #
687 # Compute characteristic polynomials and characteristic constants of a Chebyshev filter.
688 # @param P = Parameter class instance containing filter, synthesis and sweep parameters.
689 # @param FT = Frequency transformation class instance.
690 # @param symmetrizeZeros = This parameter is ignored in this class of filter.
691 # @return Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$ over characteristic constant \f$\epsilon\f$ (\f$ P(s)/\epsilon \f$).
692 # @return Es = Characteristic polynomial of the denominator of \f$S_{21}\f$ and \f$S_{11}\f$.
693 # @return Fs = Caracteristic polynomial of the numerator of \f$S_{11}\f$ considering a lossy system.
694 # @return eps = Constant in the denominator of \f$S_{21}\f$.
695 # @return epsR = Constant in the denominator of \f$S_{11}\f$.
696 # @return rootsErr = Maximum distance between the theoretical position of SQF roots and the actual position. Equal to zero for this kind of filter.
697 #
698 # \section chebyCharPolys Algorithm
699 # [TN 102.1 sec. 7.3]
700 # <ol>
701 # <li> \f$ \epsilon = \frac{1} {\sqrt{10^{L_R/10} - 1} } \f$,
702 # where \f$ L_R \f$ are the return losses (parameter "chebyLr") [TN102.1 eq.(25)].
703 #
704 # <li> \f$ P(s) = 1 \f$
705 #
706 # <li> Since \f$ E(s) \f$ has been obtained from its roots, the highest degree coeffcicient is equal to 1.
707 # In order to guarantee that
708 # \f[ \max(S_{21}(\omega)) = \max \left(\frac{P(s)}{\epsilon E(s)} \right) = 1 \f]
709 # characteristic constant \f$ \epsilon \f$ is recomputed with CharPolys.findEps() method.
710 #
711 # <li> \f$ F(s) = T_N(-js) \f$ where \f$ T_N \f$ is the Chebyshev polynomial of order N [TN 102.1 eq. (32)].
712 #
713 # The polyCheby() method is used to compute \f$ T_N \f$,
714 # while polyOmega() makes the change of variable \f$ \omega = -js \f$.
715 #
716 # <li> The roots of \f$ E(s) \f$ are [TN 102.1 eq. (34)]:
717 # \f{eqnarray*}{
718 # s_k &=& -\sinh(\mu) \sin(\theta) + j \cosh(\mu) \cos(\theta) \\
719 # \mu &=& \frac{1}{N} \sinh^{-1} \left( \frac{1}{\epsilon} \right) \\
720 # \theta &=& \left( \frac{\pi}{2} \right) \left( \frac{2k - 1}{N} \right) \\
721 # k &=& 1 \ldots N
722 # \f}
723 # The coeffcicients of \f$ E(s) \f$ are obtained from the roots using the numpy.poly1d methods.
724 #
725 # It is not necessary to define \f$ k = 1 \ldots 2N \f$ and pick up the roots in the LHS of the s-plane,
726 # as suggested in [TN 102.1 eq. (34)], since \f$ \epsilon>0 \Rightarrow \mu>0 \f$
727 # and \f$ k = 1 \ldots N \Rightarrow \sinh(\mu) \sin(\theta) > 0 \f$.
728 #
729 # <li> \f$ \epsilon_R \f$ is the ratio of the highest degree coefficients of \f$ F(s) \f$ and \f$ E(s) \f$
730 # in order to guarantee that \f$ |S_{11}(\omega)| \rightarrow 1\f$ as\f$ \omega \rightarrow \infty \f$.
731 #
732 # </ol>
733 #
734 def chebyshevFilter(P, FT, symmetrizeZeros = None):
735
736 N = P.N
737
738 # Constant eps [TN 102.1, eq. (25)]
739 eps = 1.0 / sqrt(10**(P.chebyLr/10.0) - 1)
740
741 # E(s), TN 102.1, eq. (34)
742 k = np.arange(1, N+1) # If eps>0, then mu>0 and k = 1...N yields: sinh(mu)*sin(theta) > 0
743 theta = (pi/2.0) * (2*k - 1.0)/N
744 mu = (1.0/N) * np.arcsinh(1/eps)
745 rootsEs = -sinh(mu)*np.sin(theta) + 1j*cosh(mu)*np.cos(theta)
746 Es = np.poly1d(rootsEs, r=1, variable='s')
747 Es = Es/Es.coef[0] # Normalize Es
748
749 # P(s), TN 102.1 sec. 7.3
750 Ps = np.poly1d(1, variable='s' )
751 # Orthogonailty law: [David's PFC pp. 11-12], [Cameron pp. 213]
752 if (N - Ps.o)% 2 == 0: Ps *= 1j
753
754 # Recomputation of eps to guarantee max(S21)=1
755 eps = findEps(Ps, Es)
756
757 # F(s), TN 102.1 eq. (32). Computation using recursion equation of Chebychev polynomials.
758 Tn = polyCheby(N)[N] # Chebyshev polynomial of order N.
759 Fs = polyOmega(Tn) # Change of variable w=-js
760 Fs = Fs / Fs.coef[0] # Normalize to highest degree coef. equal to 1
761
762 # epsR, TN 102.1, sec 7.3
763 epsR = Fs.coef[0] / Es.coef[0]
764 # epsR = findEps(Fs, Es)
765
766 rootsErr = 0
767
768 return Ps, Es, Fs, eps, epsR , rootsErr
769
770
771 ##
772 #
773 # Compute characteristic polynomials and characteristic constants of a Quasieliptic filter.
774 # @param P = Parameter class instance containing filter, synthesis and sweep parameters.
775 # @param FT = Frequency transformation class instance.
776 # @param symmetrizeZeros = This parameter is ignored in this class of filter.
777 # @return Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$ over characteristic constant \f$\epsilon\f$ (\f$ P(s)/\epsilon \f$).
778 # @return Es = Characteristic polynomial of the denominator of \f$S_{21}\f$ and \f$S_{11}\f$.
779 # @return Fs = Caracteristic polynomial of the numerator of \f$S_{11}\f$ considering a lossy system.
780 # @return eps = Constant in the denominator of \f$S_{21}\f$.
781 # @return epsR = Constant in the denominator of \f$S_{11}\f$.
782 # @return rootsErr = Maximum distance between the theoretical position of SQF roots and the actual position.
783 # Real or imaginary parts with abs() smaller than tolerance will be considered as a zero.
784 # A pair of roots separated less than this tolerance in s-plane will be cosidered as a double root.
785 #
786 # \section quasiElipticCharPolys Algorithm
787 # [TN 102.1 sec. 7.4]
788 # <ol>
789 # <li> \f$ \epsilon = \frac{1} {\sqrt{10^{L_R/10} - 1} } \f$,
790 # where \f$ L_R \f$ are the return losses (parameter "chebyLr") [TN102.1 eq.(25)].
791 #
792 # <li> \f$ P(s) = s^2 + a^2 \f$, [TN 102.1, eq. (40)]
793 #
794 # <li> In order to obtain \f$ E(s) \f$, the first step is to compute the following linear combination of Chebyshev polynomials [Tn 102.1, eq. (37) and (38)]:
795 # \f{eqnarray*}{
796 # \alpha &=& \frac{(a+\sqrt{a^2-1})^2}{2} \\
797 # \beta &=& -a^2 \\
798 # \gamma &=& \frac{(a-\sqrt{a^2-1})^2}{2} \\
799 # K_N(\omega) &=& \alpha \, \omega \, T_{N-1}(\omega) + \beta \, T_{N-2}(\omega) * + \gamma \, \omega \, T_{N-3}(\omega)
800 # \f}
801 # Then, \f$ E(j\omega)E(-j\omega) \f$ is [TN 102.1, eq. (41)]:
802 # \f[
803 # E(\omega)E(-\omega) = \left( a^2 - \omega^2 \right)^2 + \epsilon^2 \, K^2_N (\omega)
804 # \f]
805 # \f$ E(s)E(-s) \f$ is found by replacing \f$ \omega = -js \f$ above by the function polyOmega().
806 #
807 # Finally, the roots of \f$ E(s) \f$ are the roots of \f$ E(s)E(-s) \f$ that lie in the LHS of the s-plane.
808 #
809 # <li> Since \f$ E(s) \f$ has been obtained from its roots, the highest degree coeffcicient is equal to 1.
810 # According to eq. (39) and (41) in TN102.1, the polynomial \f$ E(s)E(-s)|_{s=j\omega} \f$
811 # is the denominator of \f$ |S_{21}(\omega)|^2 \f$ and therefore includes the constant \f$ \epsilon^2 \f$.
812 # Consequently, \f$ \epsilon \f$ must be set to the square root of the highest degree coeffcicient of \f$ E(s)E(-s) \f$.
813 #
814 # <li> \f$ F(s) \f$ follows from the energy conservation equation [TN 102.1, eq. (29)]:
815 # \f[
816 # \left| S_{11}(s) \right|^2 + \left| S_{21}(s) \right|^2 = 1
817 # \f]
818 # \f[
819 # \frac{\left| F(s) \right|^2}{\epsilon_R^2} = \left| E(s) \right|^2 - \frac{\left| P(s) \right|^2}{\epsilon^2}
820 # \f]
821 # \f[
822 # \frac{F(s)F(-s)}{\epsilon_R^2} = E(s)E(-s) - \frac{P(s)P(-s)}{\epsilon^2}
823 # \f]
824 # where we assume \f$ \epsilon_R = 1 \f$.
825 #
826 # The roots of \f$ F(s) \f$ are obtained from the roots of F(s)F(-s) by the CharPolys.factorization() method.
827 #
828 # <li> \f$ \epsilon_R \f$ must be recomputed as the ratio of the highest degree coefficients of \f$ F(s) \f$ and \f$ E(s) \f$
829 # in order to guarantee that \f$ |S_{11}(\omega)| \rightarrow 1\f$ as\f$ \omega \rightarrow \infty \f$.
830 #
831 # </ol>
832 #
833 def quasiElipticFilter(P, FT, symmetrizeZeros = None):
834
835 N = P.N
836
837 if N < 3: raise synthError, u'Quasieliptic filter is not possible for order N=2.'
838
839 # Transmission zero position.
840 a = FT.normFreq(np.array(P.qeZero))
841
842 # Check value of normalized transmission zero a>1
843 if abs(a) < 1: raise synthError, u'Quasieliptic transmission zero in the bandpass.'
844 if a < 0: a = -a
845
846 # Constant eps [TN 102.1, eq. (25)]
847 eps = 1.0 / sqrt(10**(P.qeLr/10.0) - 1)
848
849 # P(s), TN 102.1, eq. (40)
850 Ps = np.poly1d([ 1, 0, a**2], variable='s')
851 # Orthogonailty law: [David's PFC pp. 11-12], [Cameron pp. 213]
852 if (N - Ps.o)% 2 == 0: Ps *= 1j
853
854 # Coef. of linear combination of Chebyshev polynomials [Tn 102.1, eq. (38)
855 alpha = ((a+sqrt(a**2-1))**2)/2
856 beta = -a**2
857 gamma = ((a-sqrt(a**2-1))**2)/2
858
859 # Kn(w), TN 102.1, eq. (37)
860 w = np.poly1d([1, 0])
861 Tn = polyCheby(N)
862 Knw =w * Tn[N-1] * alpha + Tn[N-2] * beta + w * Tn[N-3] * gamma
863
864 # E(s)E(s), TN 102.1, eq. (41)
865 SQEw = np.poly1d([ 1, 0, -2*a**2, 0, a**4]) + eps**2 * (Knw*Knw) # E(jw)E(-jw)
866 SQEs = polyOmega(SQEw) # E(s)E(-s): Change of variable w = -js
867 rootsSQEs = SQEs.roots
868 Es = np.poly1d( rootsSQEs[np.real(rootsSQEs) < 0] , r=1, variable='s') # E(s) from roots of SQEs in left s-plane
869
870 # SQEs is unnormalized, and contains the correct constant to yield the specified return losses.
871 # Since Es has been obtained from its roots, we have lost the constant. We have to recover it and asign it to eps:
872 Es = Es/Es.coef[0]
873 eps = findEps(Ps, Es)
874
875 # F(s)F(-s), TN 102.1, eq. (29)
876 SQFs = Es * polyConj(Es) - Ps * polyConj(Ps)/eps**2
877 (Fs, F22s, rootsErr) = factorization(SQFs, 3)
878 # Fs and F22s are normalized in factorization() to highest degree coef. equal to 1
879
880 # epsR, TN 102.1, sec 7.3
881 epsR = Fs.coef[0] / Es.coef[0]
882 # epsR = findEps(Fs, Es)
883
884 return Ps, Es, Fs, eps, epsR, rootsErr
885
886 # End filter functions
887 #!!!!!!!!!!!!!!!!!!!!!!!!
888
889 #!!!!!!!!!!!!!!!!!!!!!!!!!!!
890 # Import libextrafilters
891 # Must be AT THE END to allow CIRCULAR IMPORT
892
893 # Try to import the libextrafilters. If it is not available, set importedExtraFilters flag to False.
894 try:
895 if os.path.exists('libextrafilters.pyc'): existingExtraFilters = True
896 else: existingExtraFilters = False
897 from libextrafilters import importedCheckLicense, existingCheckLicense, generalizedChebyshevFilter
898 except ImportError:
899 ##
900 # Variable set to True if the non-free libextrafilters is available. Otherwise set to False.
901 importedExtraFilters = False
902 else:
903 importedExtraFilters = True
904
905 #!!!!!!!!!!!!!!!!!!!!!!!!!!!
906 # Import mwfiltersgui
907 # Must be AT THE END to allow CIRCULAR IMPORT
908 from mwfiltersgui import PredisZeorsDlg
909
910