robust-fourier 0.0b2

Noise- and Outlier-Robust Fourier Transform with Hermite Functions in NumPy and Numba

robust_fourier

You want to compute the Fourier transform of a signal, but your signal can be corrupted by outliers? If so, this package is for you even though you will have to say goodbye to the "fast" in Fast Fourier Transform 🏃🙅‍♀️

🏗️🚧 👷👷‍♂️👷‍♀️🏗️🚧

Currently under construction. Please check back later.

⚙️ Setup and 🪛 Development

🎁 Installation

Currently, the package is not yet available on PyPI. To install it, you can clone the repository

git clone https://github.com/MothNik/robust_fourier.git

For the following commands, a Makefile is provided to simplify the process. Its use is optional, but recommended.
From within the repositories root directory, the package can be installed for normal use

# activate your virtual environment, e.g., source venv/bin/activate
make install
# equivalent to
pip install --upgrade .

or for development (with all the development dependencies)

# activate your virtual environment, e.g., source venv/bin/activate
make install-dev
# equivalent to
pip install --upgrade .["dev"]

When working in developer mode, an environment variable has to be added to run certain scripts.

ROBFT_DEVELOPER = true

🔎 Code quality

The following checks for black, isort, pyright, mypy, pycodestyle, and ruff - that are also part of the CI pipeline - can be run with

make black-check
make isort-check
make pyright-check
make mypy-check
make pycodestyle-check
make ruff-check

# or for all at once
make check

# equivalent to
black --check --diff --color ./auxiliary_scripts ./examples ./src ./tests
isort --check --diff --color ./auxiliary_scripts ./examples ./src ./tests
pyright ./auxiliary_scripts ./examples ./src ./tests
mypy ./auxiliary_scripts ./examples ./src ./tests
ruff check ./auxiliary_scripts ./examples ./src ./tests
pycodestyle ./auxiliary_scripts ./examples ./src ./tests --max-line-length=88 --ignore=E203,W503,E704

✅❌ Tests

To run the tests - almost like in the CI pipeline - you can use

make test-xmlcov  # for an XML report
make test-htmlcov  # for an HTML report

# equivalent to
pytest --cov=robust_fourier ./tests -n="auto" --cov-report=xml -x --no-jit
pytest --cov=robust_fourier ./tests -n="auto" --cov-report=html -x --no-jit

for parallelized testing whose coverage report will be stored in the file ./coverage.xml or in the folder ./htmlcov, respectively.

〰️ Hermite functions

Being the eigenfunctions of the Fourier transform, Hermite functions are excellent candidates for the basis functions for a Least Squares Regression approach to the Fourier transform. However, their evaluation can be a bit tricky.

The module hermite_functions offers a numerically stable way to evaluate Hermite functions or arbitrary order $n$ and argument - that can be scaled with a factor $\alpha$ and shifted by a constant $\mu$:

After a slight modification of the definitions in [1], the Hermite functions can be written as

with the Hermite polynomials

By making use of logarithm tricks, the evaluation that might involve infinitely high polynomial values and at the same time infinitely small Gaussians - that are on top of that scaled by an infinitely high factorial - can be computed safely and yield accurate results.

For doing so, the relation between the dilated and the non-dilated Hermite functions

and the recurrence relation for the Hermite functions

are used, but not directly. Instead, the latest evaluated Hermite function is kept at a value of either -1, 0, or +1 during the recursion and the logarithm of a correction factor is tracked and applied when the respective Hermite function is finally evaluated and stored. This approach is based on [2].

The implementation is tested against a symbolic evaluation with sympy that uses 200 digits of precision and it can be shown that even orders as high as 2,000 can still be computed even though neither the polynomial, the Gaussian nor the factorial can be evaluated for this anymore. The factorial for example would already have overflown for orders of 170 in float64-precision.

As a sanity check, their orthogonality is part of the tests together with a test for the fact that the absolute values of the Hermite functions for real input cannot exceed the value $\frac{1}{\sqrt[4]{\pi\cdot\alpha^{2}}}$.

On top of that robust_fourier comes with utility functions to approximate some special points of the Hermite functions, namely the x-positions of their

• largest root (= outermost zero),
• largest extrema in the outermost oscillation, and
• the point where they numerically fade to zero.

from robust_fourier import hermite_approx

N = 25
ALPHA = 20.0
MU = 150.0

# 1) the x-positions at which the outermost oscillation fades below machine
# precision
x_fadeout = hermite_approx.x_fadeout(
    n=ORDER,
    alpha=ALPHA,
    x_center=MU,
)
# 2) the x-positions of the largest zeros
x_largest_zero = hermite_approx.x_largest_zeros(
    n=ORDER,
    alpha=ALPHA,
    x_center=MU,
)
# 3) the x-positions of the largest extrema
x_largest_extremum = hermite_approx.x_largest_extrema(
    n=ORDER,
    alpha=ALPHA,
    x_center=MU,
)

# 4) the Gaussian approximation of the outermost oscillation ...
left_gaussian, right_gaussian = hermite_approx.get_tail_gauss_fit(
    n=ORDER,
    alpha=ALPHA,
    x_center=MU,
)
# ... which is solved for the 50% level
x_left_fifty_percent = left_gaussian.solve_for_y_fraction(y_fraction=0.5)
x_right_fifty_percent = right_gaussian.solve_for_y_fraction(y_fraction=0.5)

# 5) the Gaussian approximation is also solved for the 1% interval as a more
# realistic (less conservative) approximation of the fadeout point
x_one_percent = hermite_approx.x_tail_drop_to_fraction(
    n=ORDER,
    y_fraction=0.01,
    alpha=ALPHA,
    x_center=MU,
).ravel()

🧮 Chebyshev Polynomials

Even though the Hermite functions have some nice properties, they are not necessarily the best choice for the Fourier transform. Choosing their scaling parameter $\alpha$ can be a bit tricky. Therefore [3] suggests using Chebyshev polynomials instead. They are only defined on the interval $[-1, 1]$ and can be scaled and shifted to fit the interval $[\mu - \alpha, \mu + \alpha]$ like

for the first kind and

for the second kind. In [3] the second kind $U$ is used, but the first kind $T$ is also implemented in robust_fourier.

📖 References

• [1] Dobróka M., Szegedi H., and Vass P., Inversion-Based Fourier Transform as a New Tool for Noise Rejection, Fourier Transforms - High-tech Application and Current Trends (2017), DOI: http://dx.doi.org/10.5772/66338
• [2] Bunck B. F., A fast algorithm for evaluation of normalized Hermite functions, BIT Numer Math (2009), 49, pp. 281–295, DOI: https://doi.org/10.1007/s10543-009-0216-1
• [3] Al Marashly, O., Dobróka, M., Chebyshev polynomial-based Fourier transformation and its use in low pass filter of gravity data, Acta Geod Geophys (2024), 59, pp. 159–181 DOI: https://doi.org/10.1007/s40328-024-00444-z
• [4] Hrycak T., Schmutzhard S., Accurate evaluation of Chebyshev polynomials in floating-point arithmetic, BIT Numer Math (2019), 59, pp. 403–416, DOI: https://doi.org/10.1007/s10543-018-0738-5