libfreefilters.py

Go to the documentation of this file.

#/usr/bin/env python
# -*- coding: utf-8 -*-

## @file

##
# @package libfreefilters Free Filters Library to compute Butterworth, Chebyshev and Quasieliptic
# characteristic polynomials with minimum insertion loss. 
# @author J.M. Rius, J. Mateu, J.M. Tamayo, C. Collado, A. Padilla and J. O'Callaghan
# <p>Dpt. Signal Tehory and Communications, Universitat Politècnica de Catalunya (UPC), 
# @version 2.0.3
# @date June 11, 2013
# <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. 
# Further features were developed under contract UPC-C7767 with Thales Alenia Space España (TAS-E).  
# Contributions to the definition of the software functionality and testing have been made by Christoph Ernst, Mónica Martínez Mendoza and other 
# ESA-ESTEC personnel, and Santiago Sobrino and Luis Roglá from TAS-E.
# <p>Copyright: &copy; Universitat Politècnica de Catalunya (UPC)  2009. <br>Contact: lossyfilters@tsc.upc.edu
# <p>This program is free software; you can redistribute it and/or modify it under the terms of
# the GNU General Public License as published by the Free Software Foundation; either 
# version 3 of the License, or (at your option) any later version.
# <p>This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; 
# without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU General Public License for more details.
# <p>You should have received a copy of the GNU General Public License along with this program 
# in file "LICENSE.GPL3";  if not, download it from http://www.gnu.org/licenses/gpl-3.0.html  .
# <p>Additional permission under GNU GPL version 3 section 7:
# <p>If you modify this Program, or any covered work, by linking or combining it with the "Extra Filters Library" (libextrafilters), 
# the "Lossy Filters Library" (liblossyfilters) or the "License Check Library" (libchecklicense), or modified versions of that libraries, 
# containing parts covered by the terms described in files LICENSE.LIBEXTRAFILTERS and LICENSE.LIBLOSSYFILTERS, 
# the licensors of this Program grant you additional permission to convey the resulting work. 
# 

from math import log10, pi, sin, cos, sqrt, exp, sinh, cosh
import types
import scipy as sc
import numpy as np
import os.path
import platform
from PyQt4.QtGui import QMessageBox

from libcommonfunc import pkgNameVersion, complexStr,  myPrint,  printHeader,  synthError, askQuestion, readVerbose, warningMsg

printHeader('Successful import:' , __doc__,  importModule = True)

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
# Import libextrafilters AT THE END of the file, to avoid CIRCULAR IMPORT #
# Import mwfiltersgui AT THE END of the file, to avoid CIRCULAR IMPORT   #
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#

##
# Module global variable equal either to None or to the GUI QMainWindow class instance, containing the application main window.
# It is equal to None (default) when the calling program is the console interface (mwfiltersgui.py --nogui).
# When creating the main window, the GUI executes the libfreefilters.setGlobalVariablesLibFreeFilters() function in order to set this variable.
# The GUI methods do not use this variable. They access the mainWindow class through ParamEditDlg.mainWindow and MainWindow 'self' attributes.
mainWindow = None

##
# 
#     Function to allow the GUI set values to the global variables of libcomonfunc module.
#     This function is used in order to avoid circular imports.
#     @param variablesDict: Variable keyword argument list. Python creates a dictionary.
#     
def setGlobalVariablesLibFreeFilters(**variablesDict):    
    for key, value in variablesDict.items(): 
        globals()[key] = value

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#
# Begin class: CharPolys #

##
# 
#     Characteristic polynomials class.
#     
#     Characteristic polynomials and constants are defined as [TN 101 eq. (1)] :
#     \f[
#         S_{11}(s) = \left[ \frac{F(s)}{\epsilon_r E(s)}  \right],
#         S_{12}(s) = S_{21} = \left[ \frac{P(s)}{\epsilon E(s)}  \right],
#         S_{22}(s) = \left[ \frac{(-1)^N F(s)^*}{\epsilon_r E(s)}  \right] 
#     \f]
#     where \f$ F(s)^* = F^*(-s) \f$ is the 'Paraconjugate' polynomial of \f$ F(s) \f$ [Cameron pp. 208 footnote]
#     and is computed using polyConj function.
#     
class CharPolys(object):
    ##
    # 
    #         Compute characteristic polynomials and characteristic constants.
    #             @param P = Parameter class instance containing filter, synthesis and sweep parameters.
    #             @param  FT = Frequency transformation class instance.
    #             @param symmetrizeZeros = Automatically symmetrize Generalized Chebyshev zeros in Hz in order to compute a folded matrix with uniform Q in resonators. 
    #                                                          Equivalent to answering "Yes" in the GUI when it asks the user if he wants to symmetrize Generalized Chebyshev zeros in Hz.
    #             @return CharPolys class instance.
    #         
    def __init__(self, P,  FT,  symmetrizeZeros = None):
           
        ##
        # Characteristic polynomial \f$ P(s) \f$
        self.Ps = None
        
        ##
        # Characteristic polynomial \f$ E(s) \f$
        self.Es = None
        
        ##
        # Characteristic polynomial \f$ F(s) \f$
        self.Fs = None
        
        ##
        # Characteristic polynomial \f$ F_{22}(s) \f$
        self.F22s = None
        
        ##
        # Characteristic constant \f$ \epsilon \f$ in the denominator of \f$ S_{12} \f$ and \f$ S_{21} \f$
        self.eps = None
        
        ##
        # Characteristic constant \f$ \epsilon_R \f$ in the denominator of \f$ S_{11} \f$ 
        self.epsR = None
        
        ##
        # Roundoff error estimation in the computation of lossless characteristic polynomials, 
        # given by the rootsErr return value in CharPolys.factorization method.
        self.rootsErr = None
        
        ##
        # Value 'a' of the (s-a)/(s+a) zero/pole added in lossy filters synthesis case 3.
        self.nuqDelta = None
        
        # Filter order
        N = P.N
        
        if N < 2: raise synthError,  "Filter order must be equal or greater than 2"
        
        # Set flags for computation of Uniform Q in special topologies of TN 102.3 (cases 1, 2 or 3 )
        # Uniform Q optimization in cases 1, 2 or 3 can be performed only if TZs are symmetric in normalized frequency
        self.flagExtraFilter = self.flagSymTZ = self.flagSymTZnorm  = self.flagSymTZunorm = self.flagComplexTZ = self.flagCase1 = self.flagCase2 = self.flagCase3 = False
        
        # Standard filters ( u'Butterworth', u'Chebyshev', u'Quasieliptic' ) have no TZs or symmetric TZs  in normalized frequency
        if P.filterTransferFunc.title() in ( u'Butterworth', u'Chebyshev', u'Quasieliptic' ):
            self.flagSymTZ = True
        else:   # Only for filters in the 'Extra Filter Library', we have to find out if TZs are symmetric or not
            self.flagExtraFilter = True
            
            # Be sure that lists of zeros are really lists
            if type(P.transImagZeros) != types.ListType:  P.transImagZeros =  [ P.transImagZeros ]
            if type(P.transComplexZeros)  != types.ListType:  P.transComplexZeros =   [ P.transComplexZeros ]
            if type(P.LPcomplexZeros) != types.ListType:  P.LPcomplexZeros =  [ P.LPcomplexZeros ]
            if type(P.LPrealZeros) != types.ListType:  P.LPrealZeros =  [ P.LPrealZeros ]
            
            # Check that all TZs are outside the passband
            if np.any( np.abs(FT.normFreq(np.array( P.transImagZeros ))) < 1):
                raise synthError,  'There are transmission zeros in the passband.'
            
            # Check if TZs are symmetric in normalized and in unnormalized frequency axis.
            self.flagSymTZunorm  = self.checkUnormSymTZ(P.transImagZeros, FT)
            self.flagSymTZnorm = self.checkNormSymTZ(FT.normFreq(np.array(P.transImagZeros)))
            self.flagSymTZ = self.flagSymTZnorm or self.flagSymTZunorm
            
            # Set complex transmission zeros flag, must be False for 'Butterworth', 'Chebyshev', 'Quasieliptic'
            if P.genChebyLP or len(P.transComplexZeros) > 0: self.flagComplexTZ = True   
            
        
        # Set case1, case2 or case3 flag if transmission zeros are symmetric and not complex
        if self.flagSymTZ and not self.flagComplexTZ: # No transmission zeros or symmetric transmission zeros
            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):
                self.flagCase1 = True
            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):
                self.flagCase2 = True
            elif (N == 6) and (P.synthTech.title() == u'Prescribed Insertion Loss') and (P.nuqTech == u'kS21+pole/zero'):
                self.flagCase3 = True
        
        # If TZs are symmetric, but not in normalized frequency, we have to symmetrize them, but first ask the user if he wants to
        if self.flagExtraFilter and any([self.flagCase1, self.flagCase2, self.flagCase3]) and not self.flagSymTZnorm:
            if symmetrizeZeros is None: # The execution comes from the GUI and symmetrizeZeros command line flag is not set.
                reply = askQuestion(  
                """<center><b>Asymmetric transmission zeros:</b></center><p>
                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>
                Transmission zeros are symmetric with respect to center frequency in Hz, but not in normalized frequency.<p>
                <center><font color="red"><b>Do you want to move transmission zeros?</b></font></center>""",  
                QMessageBox.Yes | QMessageBox.No  )
                
                symmetrizeZeros = (reply == QMessageBox.Yes)
            
            if symmetrizeZeros:
                P.transImagZeros = self.symmetrizeTZ(P.transImagZeros, FT)
            else:   # The user does not want to symmetrize zeros, so we have to disable Uniform Q optimization
                self.flagCase1 = self.flagCase2 = self.flagCase3 = False
            
        # Set the roundoff error in the computation with the current precision
        self.rootsErr = sc.arcsin(sin(pi))
        
        if P.filterTransferFunc.title() ==  u'Butterworth':                         Ps, Es, Fs, eps, epsR,  rootsErr = butterworthFilter(P,  FT,  symmetrizeZeros)
        elif P.filterTransferFunc.title() ==  u'Chebyshev':                         Ps, Es, Fs, eps, epsR,  rootsErr = chebyshevFilter(P,  FT,  symmetrizeZeros)
        elif P.filterTransferFunc.title() ==  u'Quasieliptic':                    Ps, Es, Fs, eps, epsR,  rootsErr = quasiElipticFilter(P,  FT,  symmetrizeZeros)
        
        elif P.filterTransferFunc.title() ==  u'Generalized Chebyshev':    
            if importedExtraFilters is False:
                if existingExtraFilters is False:
                    raise synthError,  """Generalized Chebyshev filter is in the Library of Extra Filters,<br>
                                which is not installed in this system. 
                                <p>Contact lossyfilters@tsc.upc.edu to get the Library of Extra Filters.
                            """
                else:
                    raise synthError,  """Importing the 'Library of Extra Filters' failed. <br>
                                You have a libextrafilters.pyc file, but it cannot be imported. <br>
                                Possibly libextrafilters.pyc was compiled using a different Python version than yours.
                                <p>Your installed Python version is: %s on %s
                                <p>Read the README.html documentation or contact lossyfilters@tsc.upc.edu for assistence.
                            """ % (platform.python_version(),  platform.system())
            elif importedCheckLicense is False:
                if existingCheckLicense is False:
                    raise synthError,   """There is no check_license.pyc file in your system.
                                It should have been installed together with libextrafilters.pyc file, 
                                which actually is in your system.
                                Please follow the 'Extra Filters Library' installation instructions.
                             """
                else:
                    raise synthError,  """Importing the 'License Check Library' failed. <br>
                                You have a libchecklicense.pyc file, but it cannot be imported. <br>
                                Possibly libchecklicense.pyc was compiled using a different Python version than yours.
                                <p>Your installed Python version is: %s on %s
                                <p>Read the README.html documentation or contact lossyfilters@tsc.upc.edu for assistence.
                            """ % (platform.python_version(),  platform.system())
            else:
                Ps, Es, Fs, eps, epsR,  rootsErr = generalizedChebyshevFilter(P,  FT,  symmetrizeZeros)
            
        else: raise synthError,  u"Filter topology type '%s' not available." % P.filterTransferFunc
        
        # Update roundoff error
        self.rootsErr = max(rootsErr,  self.rootsErr)
        
        # Compute F22(s) according to [Cameron eq. (6.15)]
        F22s = (-1)**Fs.o * polyConj(Fs)
        
        # Copy local variables to class attributes
        self.Ps, self.Es, self.Fs,  self.F22s, self.eps, self.epsR = Ps, Es, Fs,  F22s, eps, epsR
    
    ##
    # 
    #         A list of transmission zeros in unnormalized frequency (Hz) is symmetrized in the normalized frequency axis.
    #         The positive zeros are averaged with the corresponding negative zeros in the normalized frequency axis.
    #         @param TZ = List of transmission zeros. Unnormalized frequency (Hz).
    #         @param FT = Frequency transform class instance.
    #         @return TZsym = List of transmission zeros, symmetrized n the normalized frequency axis. Unnormalized frequency (Hz).
    #         
    def symmetrizeTZ(self, TZ, FT):
        tol = 0
        TZ = np.array(TZ)
        TZnorm = np.array(FT.normFreq(TZ))
        TZpos = np.sort(TZnorm[np.nonzero(TZnorm>tol)]) 
        TZneg = np.sort(-TZnorm[np.nonzero(TZnorm<-tol)])
        TZnull = len(np.nonzero(np.abs(TZnorm)<=tol)[0])
        TZsymNorm = (TZpos+TZneg)/2.
        TZsymNorm = np.concatenate((np.zeros(TZnull), TZsymNorm, -TZsymNorm))
        TZsym = [ f for f in FT.unormFreq(TZsymNorm) ]
        return TZsym
    
    ##
    # 
    #         Check if list of transmission zeros (in Hz) is symmetric with respect to center frequency f0.
    #         @param TZ = List of transmission zeros in Hz.
    #         @param FT = Frequency transform class instance. It will be used to get f0.
    #         @return symmetricFlag = True is zeros are symmetric with relative error < 1e-3, False otherwise.
    #         
    def checkUnormSymTZ(self, TZ, FT):
        if len(TZ) == 0: return True
        tol = 0
        TZ = np.array(TZ)
        TZpos = np.sort(TZ[np.nonzero(TZ - FT.f0 > tol)])
        TZneg = np.sort(TZ[np.nonzero(TZ - FT.f0 < -tol)])
        if len(TZpos)==len(TZneg):
            if np.max(np.abs((TZpos+TZneg)/(2*FT.f0)-1)) < 1e-3:
                return True
        return False
    

    ##
    # 
    #             Check if list of normalized transmission zeros is symmetric with respect to zero.
    #             @param TZ = List of normalized transmission zeros.
    #             @return symmetricFlag = True is zeros are symmetric with relative error < 1e-3, False otherwise.
    #             
    def checkNormSymTZ(self, TZ):
            if len(TZ) == 0: return True
            tol = 0
            TZ = np.array(TZ)
            TZpos =np.sort(TZ[np.nonzero(TZ > tol)])
            TZneg = np.sort(TZ[np.nonzero(TZ < -tol)])
            if len(TZpos)==len(TZneg) and np.max(np.abs(TZpos+TZneg)) < 1e-3:
                return True
            else:
                return False
    

    ##
    # 
    #         Return string with polynomial coefficients separated by commas. 
    #         Rounds polynomial coefficients to given precision using complexStr() function.
    #             @param pol = Polynomial (numpy.poly1d object).
    #             @param prec = Number of significant digits (int).
    #             @return stcoef = String representation of pol.
    #         
    def polyStr(self, pol,  prec):
        stcoef = ''
        for k in range(pol.order+1):
            #re,  im = real(pol.coef[k]),  imag(pol.coef[k])
            #stre = '%.*g' % (prec,  re)  if abs(re) > 10**-prec else '0'
            #stim = '% +.*gj' % (prec,  im) if abs(im) > 10**-prec else ''
            #stcoef += stre + stim + ', '
            stcoef += complexStr(pol.coef[k],  prec) + ', '
        return stcoef[0:-2]    
    
    
    ##
    # 
    #         Save Characteristic Polynomials and constants in output file.
    #         Rounds polynomial coefficients to given precision using complexStr() function.
    #         Polynomials are of type numpy.poly1d.
    #             @param P = Parameter class instance.
    #             @param prec = Number of correct significant digits to print.
    #         
    def saveCharPolys(self, P,  prec):
        fileName = P.outDirName + '.pol'
        fout = open(fileName, 'w')
        fout.write('# %s; version %s\n' % pkgNameVersion(__doc__) )
        fout.write('# Characteristic polynomials file.\n')
        fout.write('# Polynomial coefficients are ordered from highest to lowest degree.\n\n')
        fout.write( 'F(s) = %s\n' % self.polyStr(self.Fs, prec))
        fout.write( 'P(s) = %s\n' % self.polyStr(self.Ps, prec))
        fout.write( 'E(s) = %s\n' % self.polyStr(self.Es, prec))
        
        fout.write('\n# Constant in the denominator of S12 and S21\n')
        fout.write('eps = %s\n' % complexStr(self.eps,  prec) )
        
        fout.write('\n# Constant in the denominator of S11\n')
        fout.write('epsR = %s\n' % complexStr(self.epsR,  prec) )
        fout.close()
        myPrint("Characteristic Polynomials saved to file: %s" % fileName )

# End class: CharPolys
#!!!!!!!!!!!!!!!!!!!!!!!!!!!
           
           
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Begin auxiliary functions
    
##
# 
#         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$.
#         If \f$ N \f$ is the polynomial order, \f$P_w(\omega)\f$ and \f$ P_s(s) \f$ can be expressed as:
#         \f[
#             P_w(w) = \sum_{n=0}^{N}c_{n}\omega^{n}
#         \f]
#         \f[
#             P_s(s) = [P_w(\omega)]|_{\omega=-js} = \sum_{n=0}^{N}(-j)^{n}c_{n}s^{n}
#         \f]
#         @param Pw = Input polynomial \f$ P_w(\omega) \f$.
#         @return Ps = \f$ P_s(s) \f$ .
#     
def polyOmega(Pw):
    N   = Pw.o                                                                         # Order of the polynomial
    Ps = np.poly1d(Pw.coef * pow(-1j, np.arange(N+1.)[::-1]))        # Computation of the conjugate
    return Ps             
    
        
##
# 
#     Compute Chebyshev polynomials from order 0 to N.
#     @param N = Maximum order of polynomials (int)
#     @return Tn = List containing Chebyshev polynomials. List elements of type numpy.poly1d.
#     
#     \section alfPolyCheby Algorithm
#     
#     The recursion formula for Chebyshev polynomials
#     \f[
#         T_{n}(\omega) = 2 \omega T_{n-1}(\omega) - T_{n-2}(\omega)
#     \f]
#     is used to find \f$ T_N \f$ after setting \f$ T_0 = 1 \f$ and \f$ T_1 = \omega \f$.
#     
def polyCheby(N):
    Tn = [ np.poly1d([1]),  np.poly1d([1,  0]) ]
    for k in range(2,  N+1): Tn.append( 2*Tn[1]*Tn[k-1] - Tn[k-2] )
    return Tn


##
# 
#         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$.
#         If \f$ N \f$ is the polynomial order, \f$P\f$ can be expressed as:
#         \f[
#             P(s) = \sum_{n=0}^{N} c_{n}s^{n}
#         \f]
#         We define the paraconjugate \f$P_c(s)\f$ as [Cameron pp. 208, footnote].:
#         \f[
#             P_c(s) = P(s)^* = P^*(s^*) = P^{*}(-s) = \sum_{n=0}^{N}(-1)^{n}c_{n}^*s^{n}
#         \f]
#         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
#         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$ .
#         @param P = Input polynomial \f$ P(s) \f$.
#         @return Pc = Paraconjugate polynomial \f$P_c(s)\f$.
#     
def polyConj(P):
    N   = P.o                                                                       # Order of the polynomial
    Pc = np.poly1d(np.conj(P)*pow(-1, np.arange(N+1.)[::-1]))        # Computation of the conjugate
    return Pc
    

##
# 
#         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$.
#         Polynomials are of type numpy.poly1d.
#         @param NUM = Numerator polynomial of the function to maximize.
#         @param DEN = Denominator polynomial of the function to maximize.
#         @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.
#         To do the maximization we procede as follows:
#         \f$ SQNUM(s) =|NUM(s)|^{2} = NUM(s)NUM^{*}(-s) \f$ for \f$ s=j\omega \f$.
#         \f$ SQDEN(s)=|DEN(s)|^{2} = DEN(s)DEN^{*}(-s) \f$ for \f$ s=j\omega \f$.
#         The function f to maximize is:
#         \f[ f(s)=\frac{SQNUM(s)}{SQDEN(s)}=\frac{|NUM(s)|^{2}}{|DEN(s)|^{2}} \f]
#         \f[ f'(s)=\frac{SQNUM'(s)SQDEN(s)-SQNUM(s)SQDEN'(s)}{SQDEN(s)^{2}}=0 \f]
#         The maximum of f is obtained in one of the roots of the numerator of \f$ f'(s) \f$.
#     
def findEps(NUM, DEN):
    # \f$ SQNUM(s) =|NUM(s)|^{2} = NUM(s)NUM^{*}(-s) \f$ for \f$ s=j\omega \f$.
    SQNUM = NUM * polyConj(NUM)
    # \f$ SQDEN(s)=|DEN(s)|^{2} = DEN(s)DEN^{*}(-s) \f$ for \f$ s=j\omega \f$.
    SQDEN = DEN * polyConj(DEN)
    # The function f to maximize is:
    # \f[ f(s)=\frac{SQNUM(s)}{SQDEN(s)}=\frac{|NUM(s)|^{2}}{|DEN(s)|^{2}} \f]
    # \f[ f'(s)=\frac{SQNUM'(s)SQDEN(s)-SQNUM(s)SQDEN'(s)}{SQDEN(s)^{2}}=0 \f]
    # The maximum of f is obtained in one of the roots of the numerator of \f$ f'(s) \f$.
    DERnum = SQNUM.deriv() * SQDEN - SQNUM * SQDEN.deriv()
    jOm = DERnum.r
    POL=abs(NUM(1j*np.imag(jOm))/DEN(1j*np.imag(jOm)))
    eps = max(POL)
    return eps
        

##
# 
#     Report errors in the roots position for a polynomial of the form SQF(s) = F(s)F(-s) with real coefficients.
#     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.
#           
#     The error is computed as the maximum distance between the roots theoretical and actual positions.
#         @param roots_SQF = Array of roots. Type numpy.array.
#         @return maxErr = Maximum distance between the roots theoretical and actual positions.
#     
def rootErrors(roots_SQF):
    N = len(roots_SQF)
    maxErr = 0
    for k in range(N):
        #errConj     = min( abs(roots_SQF - np.conj(roots_SQF[k])) )  # Error in the position of the conjugate
        #errOp       = min( abs(roots_SQF + roots_SQF[k]) )            # Error in the position of the opposite
        #errOpConj = min( abs(roots_SQF + np.conj(roots_SQF[k])) )  # Error in the position of the conjugate opposite
        #maxErr = max([ maxErr, errConj, errOp,  errOpConj ])
        errSym = min( np.real ( roots_SQF + roots_SQF[k] ))
        maxErr = max( maxErr, errSym )
    return maxErr
    
    
##
# 
#         Factorization of polynomial SQF(s) of order 2N as a product of two polynomials, each one of order N.
#         Polynomials are of type numpy.poly1d.
#             @param SQF  = Polynomial to factorize as a product of two.
#             @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.
#                             <br>See documentation below.
#             @return Fs, F22s = The two polynomials in which the input polynomial SQF is factorized.
#             Fs and F22s are normalized to highest order coefficient equal to 1. 
#             F22s is formed from the roots discarded in Fs, and should be equal to 
#             \f$ (-1)^N F(s)^* = (-1)^N F^*(-s) \f$ [Cameron pp. 208,212], 
#             where \f$ F(s)^* \f$ is the 'Paraconjugate' polynomial of \f$ F(s) \f$ (see function CharPolys.polyConj).
#             @return rootsErr = Maximum distance between the theoretical position of SQF roots and the actual position. 
#                                             Real or imaginary parts with abs() smaller than tolerance will be considered as a zero. 
#                                             A pair of roots separated less than this tolerance in s-plane will be cosidered as a double root.
#             
#         \section fact Factorization
#         
#         The input polynomial is factorized as:
#         \f[
#             SQF(s) = F(s)F_{22}(s)
#         \f]
#         so that when s is restricted to the imaginary axis, \f$ s=j\omega \f$, it accomplishes the following:
#         \f[
#             |SQF(j\omega)| = F(j\omega)F^{*}(j\omega), 
#             |SQF(j\omega)| = F_{22}(j\omega)F_{22}^{*}(j\omega), \forall \omega \in \Re.
#         \f]
#         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.
#         This restriction is equivalent to the fact that SQF(s) can be factorized as:
#         \f[
#             SQF(s) = F(s)F(-s)
#         \f]
#         
#         There are different methods to choose the roots of F(s) from the roots of SQF(s): <br>
#         <ol>
#             <li>If 'method' argument is an integer:
#                 <ul>
#                 <li>method == 1 -> All zeros in the left s-plane (Hurwitz).
#                 <li>method == 2 -> All zeros in the right s-plane.
#                 <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.
#                 <li>method == 4 -> The complementary of 3.
#                 </ul>
#             <li>If 'method' argument is a list of Booleans:
#                 <ul>
#                     <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).
#                     <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).
#                     <li>The zeros zero_in_right_s-plane and zero_in_left_s-plane are ordered by decreasing imaginary part.          
#                 </ul>
#         </ol>    
#         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:
#             \f[
#                 F^{1}(s) = F^{2}_{22}(s), F^{2}(s) = F^{1}_{22}(s),
#                 F^{3}(s) = F^{4}_{22}(s), F^{4}(s) = F^{3}_{22}(s).
#             \f]
#     
def factorization(SQF, method):
    
    N = SQF.o                           # Degree of the polynomial SQF to factorize.
    roots_SQF = SQF.roots       # Roots of the polynomial SQF to factorize.
    
    # Assuming quad symmetry of complex roots in SQF, test roundoff error in the positiom of roots
    rootsErr = rootErrors(roots_SQF)
    
    # Separation into different kind of roots.
    roots_PRP = roots_SQF[np.nonzero( np.real(roots_SQF) > rootsErr) ]                       # Complex roots  (not purely imaginary) with positive real part.
    roots_NRP = roots_SQF[np.nonzero( np.real(roots_SQF) < -rootsErr) ]                     # Complex roots (not purely imaginary) with negative real part.
    roots_ZRP = roots_SQF[np.nonzero( np.abs(np.real(roots_SQF)) <= rootsErr) ]      # Purely imaginary roots with zero real part (they are double roots)
    roots_ZPRP = vecSort(roots_ZRP)[::2].copy()                                                    # Purely imaginary roots with zero-positive real part (simple roots)
    roots_ZNRP = vecSort(roots_ZRP)[1::2].copy()                                                  # Purely imaginary roots with zero-negative real part (simple roots)
    
    # Remove round off errors
    roots_ZPRP = 1j*np.imag(roots_ZPRP)
    roots_ZNRP = 1j*np.imag(roots_ZNRP)
    
    # Join vectors of zeros
    roots_PRP = vecSort(np.concatenate( ( roots_PRP, roots_ZPRP ),  1))
    roots_NRP = vecSort(np.concatenate( ( roots_NRP, roots_ZNRP ),  1))
    
    if readVerbose():
        myPrint('\nFactorization of F(s)F(-s)):')
        myPrint('Zeros in left-s plane =' + str( [ complexStr(r, 3) for r in  roots_NRP ] ) )
        myPrint('Zeros in right-s plane =' + str( [ complexStr(r, 3) for r in  roots_PRP ] ) )
    
    # Open GUI for selecting zeros of F(s), only if method is already a list of bools
    if mainWindow is not None and type(method).__name__ == 'list' and type(method[0]).__name__ == 'bool':
        PZdlg = PredisZeorsDlg(roots_NRP, method, mainWindow)
        if PZdlg.exec_():
            pass
        else: 
            pass
        del PZdlg
    
    # Compute the roots for F and F22 depending on the chosen method.
    if type(method).__name__ == 'int':
        if method == 1:  # F with all roots with negative real part.
            roots_F     = vecSort(roots_NRP.copy())
            roots_F22 = vecSort(roots_PRP.copy())
        elif method == 2: # F with all roots with positive real part.
            roots_F     = vecSort(roots_PRP.copy())
            roots_F22 = vecSort(roots_NRP.copy())
        elif method == 3: 
            roots_F     = vecSort(np.concatenate( (roots_PRP[::2], roots_NRP[1::2]),  1))
            roots_F22 = vecSort(np.concatenate( (roots_PRP[1::2], roots_NRP[::2]),  1))
        elif method == 4: 
            roots_F     = vecSort(np.concatenate( (roots_PRP[1::2], roots_NRP[::2]),  1))
            roots_F22 = vecSort(np.concatenate( (roots_PRP[::2], roots_NRP[1::2]),  1))
        else:
            raise synthError,  u"Factorization of F(s)F(-s): Factorization method is not an integer in the range [1,4]."
    
    elif type(method).__name__ == 'list' and type(method[0]).__name__ == 'bool':
        if len(method) != len(roots_PRP): raise synthError,  u"Length of Factorization method boolean list is different than filter order."
        aMethod = np.asarray(method)    # Convert method to numpy array, for indexing
        roots_F     = vecSort(np.concatenate( (roots_PRP[ ~aMethod], roots_NRP[   aMethod]),  1))
        roots_F22 = vecSort(np.concatenate( (roots_PRP[    aMethod], roots_NRP[~aMethod]),  1))
            
    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."
        
    if readVerbose():
        myPrint('Factorization method = ' + str(method))
        myPrint('F(s)  =' + str( [ complexStr(r, 3) for r in  roots_F ]  ) )
        myPrint('F(-s) =' + str( [ complexStr(r, 3) for r in  roots_F22 ] ) )
        myPrint('')
    
    # Compute polynomials from their roots
    Fs = np.poly1d(roots_F, r=1)
    F22s = np.poly1d(roots_F22, r=1)
    
    # Normalize polynomials to have highest order coef equal to 1
    Fs = Fs / Fs.coef[0]
    F22s =F22s / F22s.coef[0]
    
    # Return the two polynomials taking into account that the first coeffcicient of SQF could be different from one 
    return  Fs,  F22s ,  rootsErr
        

##
# 
#             Sort a vector of complex numbers so that their imaginary parts is in decreasing order.
#                 @param vec  = complex vector to sort.
#                 @return vecsort = sorted vector.
#         
def vecSort(vec):
        # Index order for the imaginary parts in decreasing order
        ind = np.argsort(np.imag(vec))[: : -1]
        # Reorder the input vector vec
        vecsort = vec[ind]
        return vecsort  
    
    
# End: auxiliary functions
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
    
    
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Begin filter functions

##
# 
#         Compute characteristic polynomials and characteristic constants of a Butterwoth filter.
#             @param P = Parameter class instance containing filter, synthesis and sweep parameters.
#             @param  FT = Frequency transformation class instance.
#             @param symmetrizeZeros = This parameter is ignored in this class of filter.  
#             @return Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$ over characteristic constant \f$\epsilon\f$ (\f$ P(s)/\epsilon \f$).
#             @return Es = Characteristic polynomial of the denominator of \f$S_{21}\f$ and \f$S_{11}\f$.
#             @return Fs = Caracteristic polynomial of the numerator of \f$S_{11}\f$ considering a lossy system. 
#             @return eps = Constant in the denominator of \f$S_{21}\f$.
#             @return epsR = Constant in the denominator of \f$S_{11}\f$. 
#             @return rootsErr = Maximum distance between the theoretical position of SQF roots and the actual position. Equal to zero for this kind of filter.
#             
#         \section butterCharPolys Algorithm
#             [TN 102.1 sec. 7.2].
#         <ol>
#         <li> The roots of \f$ E(s)E(-s) \f$ can be found as [TN 102.1 eq. (28)]:
#        \f[
#             s_k = \left\{  \begin{array}{ll}  
#                         \exp \left( \frac{j\pi(2k-1)}{2N}  \right)   & \text{$N$ odd} \\
#                         \exp \left( \frac{j\pi 2k   }{2N}  \right)   & \text{$N$ even}  
#                     \end{array} \right.
#         \f]
#         where \f$ k = 1 \ldots 2N \f$ and \f$ N \f$ is the filter order. 
#         
#         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.
#         
#         <li> \f$ P(s) = 1 \f$
#         
#         <li> \f$ F(s) = s^N \f$
#         
#         <li> \f$ \epsilon = 1 \f$ 
#         
#         <li> \f$ \epsilon_R = 1 \f$ 
#         </ol>
#     
def butterworthFilter(P,  FT,  symmetrizeZeros = None):
    
    N = P.N

    if N>=10:  warningMsg(u'<b>WARNING:</b> <p>Butterworth filter coupling matrices are not accurate for order 10 or more.')
    
    # Roots of E(S) [TN 102.1 eq. (28)]
    k = np.arange(1,  2*N+1)
    if N%2 == 0: # N even
        rootsEs = np.exp(1j * pi  * (2*k-1) / (2*N) )
    else: # N odd
        rootsEs = np.exp(1j * pi  * k / N)
    rootsEs = rootsEs[np.where(np.real(rootsEs) < 0)]    # Pick up roots in the left plane

    #Es = np.poly1d( poly(rootsEs),  variable='s' )
    Es = np.poly1d( rootsEs,  r=1,  variable='s' )
    Es = Es/Es.coef[0]      # Normalize Es
    
    # P(s): TN 102.1 sec. 7.2
    Ps = np.poly1d(1.,  variable='s' ) 
    # Orthogonailty law: [David's PFC pp. 11-12], [Cameron  pp. 213]
    if (N - Ps.o)% 2 == 0:  Ps *= 1j    
    
    # F(s): TN 102.1 sec. 7.2
    Fs = np.poly1d(np.concatenate(( np.array([1]), np.zeros(N) )),  variable='s' )
    eps = 1.
    epsR = 1.        
        
    rootsErr = 0
    
    return Ps, Es, Fs, eps, epsR ,  rootsErr
    
    
##
# 
#         Compute characteristic polynomials and characteristic constants of a Chebyshev filter.
#             @param P = Parameter class instance containing filter, synthesis and sweep parameters.
#             @param  FT = Frequency transformation class instance.
#             @param symmetrizeZeros = This parameter is ignored in this class of filter.  
#             @return Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$ over characteristic constant \f$\epsilon\f$ (\f$ P(s)/\epsilon \f$).
#             @return Es = Characteristic polynomial of the denominator of \f$S_{21}\f$ and \f$S_{11}\f$.
#             @return Fs = Caracteristic polynomial of the numerator of \f$S_{11}\f$ considering a lossy system. 
#             @return eps = Constant in the denominator of \f$S_{21}\f$.
#             @return epsR = Constant in the denominator of \f$S_{11}\f$. 
#             @return rootsErr = Maximum distance between the theoretical position of SQF roots and the actual position. Equal to zero for this kind of filter.
#             
#         \section chebyCharPolys Algorithm
#         [TN 102.1 sec. 7.3]
#         <ol>
#         <li> \f$ \epsilon =  \frac{1} {\sqrt{10^{L_R/10} - 1} } \f$, 
#         where \f$ L_R \f$ are the return losses (parameter "chebyLr") [TN102.1 eq.(25)].
#      
#         <li> \f$ P(s) = 1 \f$
#         
#         <li> Since \f$ E(s) \f$  has been obtained from its roots, the highest degree coeffcicient is equal to 1. 
#         In order to guarantee that 
#         \f[ \max(S_{21}(\omega)) = \max \left(\frac{P(s)}{\epsilon E(s)} \right) = 1 \f]
#         characteristic constant \f$ \epsilon \f$ is recomputed with CharPolys.findEps() method.
#          
#         <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)]. 
#          
#         The polyCheby() method is used to compute \f$ T_N \f$, 
#         while polyOmega() makes the change of variable \f$ \omega = -js \f$.
#          
#         <li> The roots of \f$ E(s) \f$ are [TN 102.1 eq. (34)]:
#         \f{eqnarray*}{
#             s_k &=& -\sinh(\mu) \sin(\theta)  + j \cosh(\mu) \cos(\theta) \\
#             \mu &=& \frac{1}{N} \sinh^{-1} \left( \frac{1}{\epsilon} \right) \\
#             \theta &=& \left( \frac{\pi}{2} \right) \left( \frac{2k - 1}{N} \right) \\
#             k &=& 1 \ldots N
#         \f}
#         The coeffcicients of \f$ E(s) \f$ are obtained from the roots using the numpy.poly1d methods.
#         
#         It is not necessary to define \f$ k = 1 \ldots 2N \f$ and pick up the roots in the LHS of the s-plane, 
#         as suggested in [TN 102.1 eq. (34)], since \f$ \epsilon>0 \Rightarrow \mu>0 \f$ 
#         and \f$ k = 1 \ldots N \Rightarrow  \sinh(\mu) \sin(\theta) > 0 \f$.
#      
#         <li> \f$ \epsilon_R \f$ is the ratio of the highest degree coefficients of \f$ F(s) \f$  and \f$ E(s) \f$ 
#         in order to guarantee that \f$ |S_{11}(\omega)| \rightarrow 1\f$ as\f$ \omega \rightarrow \infty \f$.
#         
#         </ol>
#     
def chebyshevFilter(P,  FT,  symmetrizeZeros = None):
    
    N = P.N
                
    # Constant eps [TN 102.1, eq. (25)]
    eps = 1.0 / sqrt(10**(P.chebyLr/10.0) - 1)
    
    # E(s), TN 102.1, eq. (34)
    k = np.arange(1, N+1)              # If eps>0, then mu>0 and k = 1...N yields: sinh(mu)*sin(theta) > 0
    theta = (pi/2.0) * (2*k - 1.0)/N
    mu = (1.0/N) * np.arcsinh(1/eps)
    rootsEs = -sinh(mu)*np.sin(theta)  + 1j*cosh(mu)*np.cos(theta)
    Es = np.poly1d(rootsEs,  r=1,  variable='s')
    Es = Es/Es.coef[0]       # Normalize Es
    
    # P(s), TN 102.1 sec. 7.3
    Ps = np.poly1d(1,  variable='s' ) 
    # Orthogonailty law: [David's PFC pp. 11-12], [Cameron  pp. 213]
    if (N - Ps.o)% 2 == 0:  Ps *= 1j  
    
    # Recomputation of eps to guarantee max(S21)=1
    eps = findEps(Ps, Es)
    
    # F(s), TN 102.1 eq. (32). Computation using recursion equation of Chebychev polynomials.
    Tn = polyCheby(N)[N]                # Chebyshev polynomial of order N.
    Fs = polyOmega(Tn)        # Change of variable w=-js
    Fs = Fs / Fs.coef[0]        # Normalize to highest degree coef. equal to 1
    
    # epsR, TN 102.1, sec 7.3
    epsR = Fs.coef[0] / Es.coef[0]
#        epsR = findEps(Fs, Es)     
        
    rootsErr = 0
    
    return Ps, Es, Fs, eps, epsR ,  rootsErr    
    
    
##
# 
#         Compute characteristic polynomials and characteristic constants of a Quasieliptic filter.
#             @param P = Parameter class instance containing filter, synthesis and sweep parameters.
#             @param  FT = Frequency transformation class instance.
#             @param symmetrizeZeros = This parameter is ignored in this class of filter.  
#             @return Ps = Characteristic polynomial of the numerator of \f$S_{21}\f$ over characteristic constant \f$\epsilon\f$ (\f$ P(s)/\epsilon \f$).
#             @return Es = Characteristic polynomial of the denominator of \f$S_{21}\f$ and \f$S_{11}\f$.
#             @return Fs = Caracteristic polynomial of the numerator of \f$S_{11}\f$ considering a lossy system. 
#             @return eps = Constant in the denominator of \f$S_{21}\f$.
#             @return epsR = Constant in the denominator of \f$S_{11}\f$. 
#             @return rootsErr = Maximum distance between the theoretical position of SQF roots and the actual position. 
#                                             Real or imaginary parts with abs() smaller than tolerance will be considered as a zero. 
#                                             A pair of roots separated less than this tolerance in s-plane will be cosidered as a double root.
#             
#         \section quasiElipticCharPolys Algorithm
#         [TN 102.1 sec. 7.4]
#         <ol>
#         <li> \f$ \epsilon =  \frac{1} {\sqrt{10^{L_R/10} - 1} } \f$, 
#         where \f$ L_R \f$ are the return losses (parameter "chebyLr") [TN102.1 eq.(25)].
#         
#         <li> \f$ P(s) = s^2 + a^2 \f$, [TN 102.1, eq. (40)]
# 
#         <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)]:
#         \f{eqnarray*}{
#             \alpha &=& \frac{(a+\sqrt{a^2-1})^2}{2} \\
#             \beta  &=& -a^2 \\
#             \gamma &=& \frac{(a-\sqrt{a^2-1})^2}{2}   \\
#             K_N(\omega) &=& \alpha \, \omega \, T_{N-1}(\omega) + \beta \, T_{N-2}(\omega) * + \gamma \, \omega \, T_{N-3}(\omega)
#         \f} 
#         Then, \f$ E(j\omega)E(-j\omega) \f$ is [TN 102.1, eq. (41)]:
#         \f[
#             E(\omega)E(-\omega) = \left( a^2 - \omega^2 \right)^2 + \epsilon^2 \, K^2_N (\omega)
#         \f]
#         \f$ E(s)E(-s) \f$ is found by replacing \f$ \omega = -js \f$ above by the function polyOmega().
#         
#         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.
#         
#         <li> Since \f$ E(s) \f$  has been obtained from its roots, the highest degree coeffcicient is equal to 1. 
#         According to eq. (39) and (41) in TN102.1, the polynomial \f$ E(s)E(-s)|_{s=j\omega} \f$ 
#         is the denominator of \f$ |S_{21}(\omega)|^2 \f$ and therefore includes the constant \f$ \epsilon^2 \f$. 
#         Consequently, \f$ \epsilon \f$ must be set to the square root of the highest degree coeffcicient of \f$ E(s)E(-s) \f$.
#         
#         <li> \f$ F(s) \f$ follows from the energy conservation equation [TN 102.1, eq. (29)]:
#             \f[
#             \left| S_{11}(s) \right|^2 + \left| S_{21}(s) \right|^2 = 1
#             \f]
#             \f[
#             \frac{\left| F(s) \right|^2}{\epsilon_R^2} = \left| E(s) \right|^2 - \frac{\left| P(s) \right|^2}{\epsilon^2}
#             \f]
#             \f[
#             \frac{F(s)F(-s)}{\epsilon_R^2} = E(s)E(-s) - \frac{P(s)P(-s)}{\epsilon^2}
#             \f]
#             where we assume \f$ \epsilon_R = 1 \f$.
#             
#             The roots of \f$ F(s) \f$ are obtained from the roots of F(s)F(-s) by the CharPolys.factorization() method. 
#             
#             <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$ 
#         in order to guarantee that \f$ |S_{11}(\omega)| \rightarrow 1\f$ as\f$ \omega \rightarrow \infty \f$.
#         
#         </ol>
#     
def quasiElipticFilter(P,  FT,  symmetrizeZeros = None):
    
    N = P.N
    
    if N < 3: raise synthError,  u'Quasieliptic filter is not possible for order N=2.'
    
    # Transmission zero position.
    a = FT.normFreq(np.array(P.qeZero))

    # Check value of normalized transmission zero a>1
    if abs(a) < 1: raise synthError,  u'Quasieliptic transmission zero in the bandpass.'
    if a < 0: a = -a
    
    # Constant eps [TN 102.1, eq. (25)]
    eps = 1.0 / sqrt(10**(P.qeLr/10.0) - 1)
    
    # P(s), TN 102.1, eq. (40)
    Ps  = np.poly1d([ 1,  0,  a**2],  variable='s')
    # Orthogonailty law: [David's PFC pp. 11-12], [Cameron  pp. 213]
    if (N - Ps.o)% 2 == 0:  Ps *= 1j  
    
    # Coef. of linear combination of Chebyshev polynomials [Tn 102.1, eq. (38)
    alpha = ((a+sqrt(a**2-1))**2)/2
    beta  = -a**2
    gamma = ((a-sqrt(a**2-1))**2)/2
    
    # Kn(w), TN 102.1, eq. (37)
    w = np.poly1d([1,  0])
    Tn = polyCheby(N)
    Knw =w * Tn[N-1] * alpha + Tn[N-2] * beta + w * Tn[N-3] * gamma
    
    # E(s)E(s), TN 102.1, eq. (41)
    SQEw = np.poly1d([ 1,  0,  -2*a**2,  0,  a**4]) + eps**2 * (Knw*Knw)         # E(jw)E(-jw)
    SQEs = polyOmega(SQEw)          # E(s)E(-s): Change of variable w = -js
    rootsSQEs = SQEs.roots
    Es = np.poly1d( rootsSQEs[np.real(rootsSQEs) < 0] ,  r=1,  variable='s')     # E(s) from roots of SQEs in left s-plane
    
    # SQEs is unnormalized, and contains the correct constant to yield the specified return losses.
    # Since Es has been obtained from its roots, we have lost the constant. We have to recover it and asign it to eps: 
    Es = Es/Es.coef[0] 
    eps = findEps(Ps, Es) 
    
    # F(s)F(-s), TN 102.1, eq. (29)
    SQFs = Es * polyConj(Es) - Ps * polyConj(Ps)/eps**2
    (Fs,  F22s,  rootsErr) = factorization(SQFs,  3)
    # Fs and F22s are normalized in factorization() to highest degree coef. equal to 1
        
    # epsR, TN 102.1, sec 7.3
    epsR = Fs.coef[0] / Es.coef[0]
#    epsR = findEps(Fs, Es)     
    
    return Ps, Es, Fs, eps, epsR,  rootsErr
    
# End filter functions
#!!!!!!!!!!!!!!!!!!!!!!!!

#!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Import libextrafilters
# Must be AT THE END to allow CIRCULAR IMPORT

# Try to import the libextrafilters. If it is not available, set importedExtraFilters flag to False.
try:
    if os.path.exists('libextrafilters.pyc'): existingExtraFilters = True
    else: existingExtraFilters = False
    from libextrafilters import  importedCheckLicense,  existingCheckLicense,  generalizedChebyshevFilter
except ImportError: 
    ##
    # Variable set to True if the non-free libextrafilters is available. Otherwise set to False.
    importedExtraFilters = False     
else: 
    importedExtraFilters = True    
    
#!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Import mwfiltersgui
# Must be AT THE END to allow CIRCULAR IMPORT    
from mwfiltersgui import PredisZeorsDlg