将一段八度（mrdivide）代码转换为numpy

Question

我试图将一段倍频程代码转换为numpy，但是我从倍频程和numpy得到了不同的结果。以下是我的数据（事实上，它们远远大于下面给出的数据）：

A =

 Columns 1 through 6:

   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.908807141567176e-01   9.954090291002812e-01   9.908807141567179e-01
   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01
   9.908807141567176e-01   9.776578008378289e-01   9.567491788333946e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01
   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01
   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00
   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.608235911808946e-01   9.820064469806473e-01   9.954090291002812e-01
   9.776578008378289e-01   9.649505047327671e-01   9.448300707857176e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01
   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01
   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01
   9.567491788333946e-01   9.448300707857176e-01   9.259179407344914e-01   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01
   9.608235911808946e-01   9.567491788333946e-01   9.448300707857176e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01
   9.567491788333946e-01   9.608235911808946e-01   9.567491788333946e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01
   9.448300707857176e-01   9.567491788333946e-01   9.608235911808946e-01   9.448300707857176e-01   9.649505047327673e-01   9.776578008378289e-01

 Columns 7 through 12:

   9.776578008378289e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01   9.567491788333946e-01
   9.908807141567179e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.448300707857176e-01
   9.954090291002812e-01   9.448300707857176e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.259179407344914e-01
   9.608235911808946e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01   9.820064469806473e-01
   9.820064469806473e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.776578008378289e-01
   9.954090291002812e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.649505047327671e-01
   1.000000000000000e+00   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.448300707857176e-01
   9.567491788333946e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.608235911808946e-01   9.954090291002812e-01
   9.776578008378289e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.908807141567179e-01
   9.908807141567179e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.776578008378289e-01
   9.954090291002812e-01   9.608235911808946e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.567491788333946e-01
   9.448300707857176e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01   1.000000000000000e+00
   9.649505047327673e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.954090291002812e-01
   9.776578008378289e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.820064469806473e-01
   9.820064469806473e-01   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.608235911808946e-01

 Columns 13 through 15:

   9.608235911808946e-01   9.567491788333946e-01   9.448300707857176e-01
   9.567491788333946e-01   9.608235911808946e-01   9.567491788333946e-01
   9.448300707857176e-01   9.567491788333946e-01   9.608235911808946e-01
   9.776578008378289e-01   9.649505047327671e-01   9.448300707857176e-01
   9.820064469806473e-01   9.776578008378289e-01   9.649505047327673e-01
   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01
   9.649505047327673e-01   9.776578008378289e-01   9.820064469806473e-01
   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01
   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01
   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01
   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.954090291002812e-01   9.820064469806473e-01   9.608235911808946e-01
   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01
   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01
   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00

以及

b =

  -1.024208397018539
  -1.055718555015945
  -1.066560607689640
  -1.187368188387253
  -1.258866007703282
  -1.305258462589997
  -1.321354530870290
  -1.333661132027602
  -1.421384660329320
  -1.478743779481671
  -1.498725636719488
  -1.385960967135295
  -1.479663779776475
  -1.541078471216082
  -1.562500000000000

在八度音阶中我有x = b'/A。我找不到/对应的numpy函数。到目前为止，我尝试了numpy的x = np.dot(b.T,np.linalg.inv(A))，但结果与八度音阶不同。你知道吗

x = b'/A的倍频程结果如下：

x =

 Columns 1 through 6:

  -5.642309525591432e+00   7.813412870218545e+00  -1.559855155426489e+02  -3.597241224212262e+01   2.201914551287831e+02  -3.100354445411479e+02

 Columns 7 through 12:

   7.253956488595386e+02   7.369595892720794e+01  -4.313273469816049e+02   6.064725968037579e+02  -9.855235323530542e+02  -4.111380598448122e+01

 Columns 13 through 15:

   2.334109900297194e+02  -3.269547704109582e+02   4.254317069619117e+02

numpy的结果是

x=np.array([[-5.642310487492068,    7.813414778371225,  -155.9855165364133, -35.9724138597885,  220.1914623928024,  -310.0354544342263, 725.3956532399461,  73.69596218669903,  -431.3273588509765, 606.472611254314,   -985.5235383154941, -41.11380770278629, 233.4109958125337,  -326.9547770833597, 425.4317096135928]])

如果有人能帮我找到和numpy的八度音阶相同的结果，我将不胜感激。或者比当前结果更接近倍频程结果。你知道吗

Answer 1

看起来Matlab和numpy的计算精度不同。我会尽量让两者都使用浮点精度。你知道吗

默认情况下，matlab使用16位精度进行数值计算（请参见https://www.mathworks.com/help/symbolic/increase-precision-of-numeric-calculations.html）。在这个网站上，它也描述了如何改变精度。你知道吗

就我而言，numpy使用32位精度（您可能会发现本文很有用：how set numpy floating point accuracy?）。你知道吗

Answer 2

获取两个数组（从文章中复制并粘贴）

x = [-5.642309525591432e+00, 7.813412870218545e+00, -1.559855155426489e+02, -3.597241224212262e+01, 2.201914551287831e+02, -3.100354445411479e+02, 7.253956488595386e+02, 7.369595892720794e+01, -4.313273469816049e+02, 6.064725968037579e+02, -9.855235323530542e+02, -4.111380598448122e+01, 2.334109900297194e+02, -3.269547704109582e+02, 4.254317069619117e+02]
y = [-5.642310487492068,    7.813414778371225,  -155.9855165364133, -35.9724138597885,  220.1914623928024,  -310.0354544342263, 725.3956532399461,  73.69596218669903,  -431.3273588509765, 606.472611254314,   -985.5235383154941, -41.11380770278629, 233.4109958125337,  -326.9547770833597, 425.4317096135928]

通过allclose运行它们

import numpy as np
np.allclose(x, y) # True

它们（在some tolerance内）是相同的。你知道吗