Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
mammalwong opened this issue Jul 16, 2024 · 1 comment

Comments

@mammalwong
Copy link

mammalwong commented Jul 16, 2024

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
Copy link
Author

mammalwong commented Jul 16, 2024

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()

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant