`iczt(X)` is numerical unstable and return nan when `len(X)` is too big (arround > 5000), fix suggestion included

[mammalwong]

version of czt: 0.0.7, issue was found when using google colab. To reproduce the issue:

import numpy as np
import czt as czt
for n in [1000,10000]:
  x = np.linspace(-4,4,n)
  y = np.exp(-x**2)
  yhat = czt.czt(y,simple=False)
  yphi = czt.iczt(yhat,simple=False)
  print(np.any(np.isnan(yhat)),np.any(np.isnan(yphi)))

outputs:

False False
False True

expected result:

False False
False False

[mammalwong]

suggested fix (tested in google colab locally) The np.cumprod call in the below lines of the source of iczt caused this numerical instability issue, it makes many tailing entries in p become 0:

    p = np.r_[1, (W ** k[1:] - 1).cumprod()]
    u = (-1) ** k * W ** (k * (k - n + 0.5) + (n / 2 - 0.5) * n) / p
    # equivalent to:
    # u = (-1) ** k * W ** ((2 * k ** 2 - (2 * n - 1) * k + n * (n - 1)) / 2) / p
    u /= p[::-1]

it can be solved by modifying the above few lines like this:

    u = (-1) ** k * W ** (k * (k - n + 0.5) + (n / 2 - 0.5) * n) # /p is removed from here (1)
    p = np.r_[1, W ** k[1:] - 1]
    lp = np.abs(p) # lp stand for ln(p)
    lp = np.cumsum(np.log(lp)) + np.angle(np.cumprod(p/lp))*1j
    # above seperate the magnitude and angle of the entries in p
    # it cumsum magnitude in log domain to replace the unstable cumprod of the magnitude
    # and cumprod only the normalized angle to ensure the angle is also stable when X is long
    u /= np.exp(lp+lp[::-1]) # /p from (1) is moved to here as +lp in exp()

[garrettj403]

Thanks for point this out and sorry it took so long to get back to you.

I added your code to commit 96fdc262238bcccd703da92fc1eb99b1f58bcd02.

[mammalwong]

Thanks for the acknowledgement.

在 2024年12月20日週五 15:13，John Garrett @.***> 寫道：

Thanks for point this out and sorry it took so long to get back to you.

I added your code to commit 96fdc26 https://github.com/garrettj403/CZT/commit/96fdc262238bcccd703da92fc1eb99b1f58bcd02 .

— Reply to this email directly, view it on GitHub https://github.com/garrettj403/CZT/issues/17#issuecomment-2556420219, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGY3BIYT757EHPLFXSIO2M32GO7TBAVCNFSM6AAAAABK66R42KVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKNJWGQZDAMRRHE . You are receiving this because you authored the thread.Message ID: @.***>