如何提高FosterCauer变换的scipy.signal.invres（）的数值稳定性？

Question

概述

我目前正在实现Foster-Cauer转换（https://arxiv.org/ftp/arxiv/papers/0801/0801.0817.pdf，第3章）。变换的一部分是对我使用^{}的部分分式分解求逆

将此函数应用于我的测试数据集（512个summands），返回r和p的以下值：

r =[1.81758254e+005 1.73495508e+010 8.18928046e+014 2.54996869e+019
 5.89531786e+023 1.07985536e+028 1.63302037e+032 2.09775523e+036
 2.33738380e+040 2.29542819e+044 2.01208717e+048 1.59048788e+052
 1.14339548e+056 7.52908473e+059 4.56887352e+063 2.56848686e+067
 1.34380079e+071 6.56947543e+074 3.01172487e+078 1.29889143e+082
 5.28501907e+085 2.03403852e+089 7.42214945e+092 2.57328050e+096
 8.49340621e+099 2.67359361e+103 8.03985241e+106 2.31316334e+110
 6.37658571e+113 1.68643162e+117 4.28434117e+120 1.04673144e+124
 2.46202142e+127 5.58077321e+130 1.22027393e+134 2.57615017e+137
 5.25538381e+140 1.03682514e+144 1.97972229e+147 3.66111575e+150
 6.56190651e+153 1.14060777e+157 1.92398223e+160 3.15123170e+163
 5.01439470e+166 7.75619953e+169 1.16679575e+173 1.70792839e+176
 2.43376135e+179 3.37766023e+182 4.56742875e+185 6.02041195e+188
 7.73843224e+191 9.70329489e+194 1.18736931e+198 1.41842941e+201
 1.65475742e+204 1.88586127e+207 2.10025981e+210 2.28643312e+213
 2.43386221e+216 2.53402551e+219 2.58121935e+222 2.57309196e+225
 2.51082238e+228 2.39892821e+231 2.24473979e+234 2.05762363e+237
 1.84806681e+240 1.62674265e+243 1.40366782e+246 1.18753464e+249
 9.85267677e+251 8.01817025e+254 6.40168198e+257 5.01525611e+260
 3.85613785e+263 2.91038543e+266 2.15656756e+269 1.56914962e+272
 1.12131239e+275 7.87082588e+277 5.42767831e+280 3.67768916e+283
 2.44888537e+286 1.60272215e+289 1.03111400e+292 6.52190767e+294
 4.05622071e+297 2.48087889e+300 1.49239357e+303 8.83101221e+305
             inf             nan             nan             nan
             ...             nan             nan             nan] 

p = [1.00000000e+000 9.96908839e+004 4.88722027e+009 1.57350549e+014
 3.74788002e+018 7.05166051e+022 1.09264854e+027 1.43512931e+031
 1.63204727e+035 1.63328269e+039 1.45699278e+043 1.17070111e+047
 8.54616211e+050 5.70928590e+054 3.51207485e+058 2.00002171e+062
 1.05929039e+066 5.23940044e+069 2.42889753e+073 1.05877375e+077
 4.35236476e+080 1.69166618e+084 6.23167492e+087 2.18040672e+091
 7.26063674e+094 2.30519691e+098 6.98983805e+101 2.02734005e+105
 5.63264124e+108 1.50108507e+112 3.84192182e+115 9.45471542e+118
 2.23964737e+122 5.11195914e+125 1.12535707e+129 2.39156913e+132
 4.91063101e+135 9.75001200e+138 1.87335208e+142 3.48574546e+145
 6.28540875e+148 1.09905071e+152 1.86474430e+155 3.07182462e+158
 4.91581394e+161 7.64630406e+164 1.15661721e+168 1.70225590e+171
 2.43872720e+174 3.40253282e+177 4.62521246e+180 6.12820557e+183
 7.91739215e+186 9.97806697e+189 1.22712280e+193 1.47320420e+196
 1.72711601e+199 1.97791985e+202 2.21342296e+205 2.42116169e+208
 2.58950884e+211 2.70876576e+214 2.77209469e+217 2.77616587e+220
 2.72143485e+223 2.61202021e+226 2.45521011e+229 2.26067725e+232
 2.03951713e+235 1.80323862e+238 1.56282881e+241 1.32798931e+244
 1.10660519e+247 9.04468128e+249 7.25239370e+252 5.70610562e+255
 4.40604729e+258 3.33954943e+261 2.48503429e+264 1.81575454e+267
 1.30297083e+270 9.18408285e+272 6.35959443e+275 4.32695990e+278
 2.89308908e+281 1.90121098e+284 1.22814874e+287 7.79982503e+289
 4.87070338e+292 2.99109181e+295 1.80657515e+298 1.07331267e+301
 6.27328293e+303 3.60757324e+306             inf             inf
             inf             ...             inf             inf
             inf]

我知道inf值是由指数增加引起的

---

问题

创建这些inf值是因为^{}背后的算法是病态的吗？ 除了^{}还有别的选择吗

测试数据

数据可在此处下载：

https://depot.tu-dortmund.de/ran43

链接有效期至：03.12.19

---

示例代码

您可以使用以下方法再现错误：

import numpy as np 
from scipy.signal import invres


def foster_cauer_transformation(R_foster, C_foster):
    r, p = invres(1/C_foster, -1/(R_foster * C_foster), k=np.zeros(1))
    # This is not the end of the foster cauer transformation, but enough for this example. 
    return p, q


# file is attached to post, see Test data section
R_foster, C_foster = np.genfromtxt('R_and_C_foster.txt', unpack=True)

r_result, p_results = foster_cauer_transformation(R_foster, C_foster)

print('r', r_result)
print('\n ------ \n')
print('p', p_result)