6
$\begingroup$

Issues with the Fourier Transform of $f(t)=(1-t^2)^4$ on $[-1,\,1]$, should be analytical but looks like having a singularity with noise-like rippling


Intro

I was trying to made a compact-supported approximation of a Gaussian Envelope $$g(t)= e^{-4t^2} \tag{Eq. 1}\label{Eq. 1}$$ So I was trying to use instead: $$f(t)= (1-t^2)^4\cdot\theta(1-t^2) \cong \left(\frac{1-t^2+|1-t^2|}{2}\right)^4 \tag{Eq. 2}\label{Eq. 2}$$ where $\theta(t)$ is the Heaviside step function.

Both which can be seen here:

gaussian envelope and compact-supported approximation

As is shown in the image, at "my taste" the match "looks good", but accurately speaking, it fulfill that:

  1. The function $f(t) = 0,\,|t|\geq 1$, so it is of compact support.
  2. At the edges of the domain $\partial t =\{-1;\,1\}$ the function $s(t)=(1-t^2)^4$ is already zero: $s(-1)=s(1)=0$, so shouldn't be much issues related with the term $\theta(1-t^2)$, since it is going to be at most an avoidable discontinuity: $\lim\limits_{t \to \partial t^{\pm}} f(t) = \lim\limits_{t \to \partial t} f(t) = 0$, so the function is continuous on $\mathbb{R}$.
  3. Using Wolfram-Alpha, the following norms shows to be bounded: $\|f\|_\infty = 1 < \infty$, $\|f\|_1 = \frac{256}{315} \approx 0.813 < \infty$, $\|f\|_2^2 = \frac{65536}{109395} \approx 0.599 < \infty$ so is bounded, absolute integrable, and energy limited.
  4. Also its derivative is: $$\begin{array}{r c l} \frac{d}{dt}\left((1-t^2)^4\cdot\theta(1-t^2)\right) & = & \frac{d}{dt}\left((1-t^2)^4\right)\cdot\theta(1-t^2)+(1-t^2)^4\cdot\frac{d}{dt}\left(\theta(1-t^2)\right)\\ & = & -8\,t\,(1-t^2)^3\cdot\theta(1-t^2)+(1-t^2)^4\cdot\frac{d}{dt}\left(\frac{1+\text{sgn}(1-t^2)}{2}\right)\\ & = & -8\,t\,(1-t^2)^3\cdot\theta(1-t^2)+2\,t\,(1-t^2)^4\cdot\delta(1-t^2)\\ & = & -8\,t\,(1-t^2)^3\cdot\theta(1-t^2)+2\,t\,(1-t^2)^3\cdot\underbrace{(1-t^2)\cdot\delta(1-t^2)}_{=\,0,\,\text{since}\,(1-x)\delta(1-x)=0} \\ \end{array}$$ $$\Rightarrow \frac{df(t)}{dt} = -8\,t\,(1-t^2)^3\cdot\theta(1-t^2) \tag{Eq. 3}\label{Eq. 3} $$ so its derivative is also continuous by the same reasons of point (ii), and $\|f'\|_\infty = \frac{1798}{343\sqrt{7}} \approx 1.9 < \infty$, $\|f'\|_1 = 2 < \infty$, and $\|f'\|_2^2 = \frac{131072}{45045} \approx 2.9 < \infty$.

So with all these point the approximation $f(t)$ look quite "well-behaved": note that is not a "smooth function", but same analysis of point $(4)$ could be done successfully for $f''$ also.

From now on I will continue labeling the observations so you can accurately point where I am having a misconception.

  1. Since the function $f(t)$ is of compact-support and also is squared-integrable (since $\|f\|_2^2 < \infty$ as show in point $(3)$), I where expecting that the Paley–Wiener theorem to be fulfilled, so the Fourier Transform of $f(t)$ was going to be an analytical function.

Using Wolfram Alpha I calculate the Fourier Transform (under the "electricians" definition), in 3 different ways showing each of them the same result so I believe is right calculated -way 1, way 2, and way 3:

$$\hat{f}(w) = \int\limits_{-\infty}^{\infty} f(t)\,e^{-iwt}\,dt = \int\limits_{-1}^{1} (1-t^2)^4\,e^{-iwt}\,dt = \frac{768 \Big(5w\,(2 w^2-21)\cos(w)+(w^4-45w^2+105) \sin(w)\Big)}{w^9} \tag{Eq. 4}\label{Eq. 4}$$

So far so good, since \eqref{Eq. 4} looks like the classical Fourier Transform made by a polynomial mixed with some trigonometric functions, but when trying to calculate $\|\hat{f}\|_1$ and $\|iw\hat{f}\|_1$ the integrals get stuck on Wolfram-Alpha, so I decide to plot it to see what is going on.

I am going to graph it against the Fourier Transform of $g(t)$, which under the same definition of the transform is going to be: $$\hat{g}(w) = \frac{\sqrt{\pi}}{2}e^{-\frac{w^2}{16}}\tag{Eq. 5}\label{Eq. 5}$$

Comparison of the Fourier Transforms

They look similar, but something "weird" is happening near the DC component $w \approx 0$ for the solution $\hat{f}(w)$ of \eqref{Eq. 4}.

To have a better insight I plot $\hat{f}(w)$ also in Desmos and there exist an "horrible noise-like rippling" that also "looks" to be diverging as a singularity!

Horrible Rippling

I tried to find the value of $\|\hat{f}\|_\infty$ with Wolfram-Aplha but I believe are numerical errors since the value change when changing the domain limits, but are numbers of the order of $10^{93}$!!!.

  1. I believe that a Analytical function is smooth (as it should be the Fourier Transform $\hat{f}$ following point $(5)$), so it should be bounded on every closed interval of its domain: maybe here I am mistaken, but since analytic functions are smooth, they should be continuous on every closed interval of the domain (since are already differentiable), and since is continuous on a closed interval it should achieve a bounded maximum and minimum because of the Extreme value theorem.

But even if I am wrong, I wasn't expecting an analytical function to be so "ill-behaved" as is being $\hat{f}$ on $w \in [-0.2,\,0.2]$.


Questions

So the main questions are:

Q1: Is $\hat{f}(w)$ of \eqref{Eq. 4} properly obtained? Also related quantities as $\|\hat{f}\|_1$ and $\|\hat{f}\|_{\infty}$

Q2: If right, Why is this "rippling" happening? It is right? or is a miscalculation of the graphic software? or a numerical issue? Actually its looks like "Brownian" or "White Noise".

Q3: Is truly having a singularity at some angular frequency $w$? (doing limits in Wolfram Alpha doesn't work, at least it is saying that $\lim_{w\to 0} \hat{f}(w)=\frac{256}{315}\approx 0.813$ which is near the main tendency, but surely below the rippling peaks - maybe the singularity is near but not in $w=0$).

Q4: Is this rippling a known result like the "Gibbs' Phenomenon? How is called? If it is a known phenomenon I would like to search for any references - maybe the is a way to avoid it.

Q5: Where I am misleading my analysis? (since this singularity shouldn't exists under my assumptions and my actual knowledge - which is quite basic by the way)


Motivation

After seeing this video about the absorption of light, temporal dispersion, and Kramers-Kronig relations, required for the complex-valued modeling of the refraction index to preserve causality (is quite a good video, simple and short), I see why dispersion must happen or other way the absorbed frequency will lead to a "time-wide" signal which will become non-causal in the analysis (it will be entering the absorptive medium before the frequency-component subtraction had happen on the first place). But since a Gaussian Envelope is used, the "causality-issue" where already present since this function only vanishes at infinity, but since spectra decays faster than $1/|w|$, the Kramers-Kronig relation still holds "hidding" this issue.

This is why I was trying to find another envelope but with an accurate compact-support, since it is multiplied by a trigonometric function it Fourier Transform is only a displaced version of the Transform of the Envelope function, and here is where I found the "problem" since the approximation I take has this "issue" in the spectrum - don't knowing now if its a "real-life issue" or just a numeric problem of the calculation algorithms.

Hope you find it has interesting I am, beforehand thanks you very much.


added later

I have noted now another interesting thing:

Integrating $\hat{f}$ as is shown here, and integrating $\hat{g}$ as is shown here gives the same result: $$\int\limits_{-\infty}^{\infty}\hat{f}(w)\,dw = \int\limits_{-\infty}^{\infty}\hat{g}(w)\,dw = 2\pi$$


I have tried another slightly different approximation: $$q(t) = \left(\frac{1-t^2+|1-t^2|}{2}\right)^{\pi} \tag{Eq. 6}\label{Eq. 6}$$

And following Wolfram-Alpha its Fourier Transform is:

$$\hat{q}(w) = \int\limits_{-\infty}^{\infty} q(t)\,e^{-iwt}\,dt = \sqrt{\pi}\,\Gamma(1+\pi)_0\tilde{F}_1\left(;\,\frac{3}{2}+\pi;\,-\frac{w^2}{4}\right) \tag{Eq. 7}\label{Eq. 7}$$ which plot doesn't show the noise-like-rippling-singularity, so with this, as it was noted on the comments and answers, the situation was a numerical issue. Even so, if you directly replace all the numbers $\pi$ by the number $4$ in \eqref{Eq. 7} the rippling appears again.

Later I noted that the result of \eqref{Eq. 4} is equivalent to: $$\hat{f}(w) = \sqrt{\pi}\,\Gamma(1+4)_0\tilde{F}_1\left(;\,\frac{3}{2}+4;\,-\frac{w^2}{4}\right) \tag{Eq. 8}\label{Eq. 8}$$

So obviously again the numerical issues arises, but for a "good" approximation I am using on Wolfram-Alpha: $$\hat{f}^*(w) = \sqrt{\pi}\,\Gamma(1+4)_0\tilde{F}_1\left(;\,\frac{3}{2}+4.0001;\,-\frac{w^2}{4}\right) \tag{Eq. 9}\label{Eq. 9}$$

and suddenly the problem is gone... it is kind of ridiculous that "just" for the exact numbers the problem happen.

Approximated version

Since the regularized confluent hypergeometric function is not quite extended on every software library I find useful the following property shown by Wolfram Alpha to make approximated plots of these kind of functions, by taking an arbitrary parameter $a$: $$\sqrt{\pi}\,\Gamma(1+a)_0\tilde{F}_1\left(;\,\frac{3}{2}+a;\,-\frac{w^2}{4}\right) = \text{sgn}(w)\sqrt{\pi}\,\left(\frac{2}{w}\right)^{a+\frac{1}{2}}\Gamma(1+a)J_{a+\frac{1}{2}}(w) \tag{Eq. 10}\label{Eq. 10}$$ with $J_n(x)$ the Bessel function of the first kind.

Unfortunately, I wasn't able to find $\|\hat{f}\|_1$ under any of these approximations (both shown on the last image), "as it where diverging" (I don't know if this is the case), even when $\|\hat{g}\|_1=2\pi << \infty$, which is counter-intuitive for me giving its similarity (the lobes should be adding even less area than a Gaussian in principle), but comparing their graphs shows that the compact-support function has lobes much more spread than the Gaussian, as is shown here.

comparison of spectrum

But nevertheless, the tails are decreasing faster than $|\frac{1}{w^2}|$ so I don´t know why the result is not finite (see here).

Tails comparison

As example, if I used the other aporoximarion: $$\hat{f}^*(w) \approx \frac{962.612\,J_{4.5001}(w)}{w^{4.5001}}\tag{Eq. 11}\label{Eq. 11}$$ Wolfram-Alpha is unable to find $\|\hat{f}^*(w)\|_1$, but it is possible to bound it as: $$\|\hat{f}^*(w)\|_1 < \underbrace{\int\limits_{-10}^{10}\left|\frac{962.612\,J_{4.5001}(w)}{w^{4.5001}}\right|\,dw}_{6.34416}+\underbrace{\int\limits_{10}^{\infty}\left|\frac{2}{|w^{2}|}\right|\,dw}_{0.2} <\infty$$ as can be seen here and here. And since the difference with the approximation and the main function with numerical issues is in the bounded domain $[-1,\,1]$, it should also be bounded, so the norm $\|\hat{f}(w)\|_1<\infty$ but somehow Wolfram-Alpha can't find it.

$\endgroup$
9
  • $\begingroup$ Are you sure about the denominator $w^9?$ Certainly, if should be $$\hat f(0)=\int_{-1}^1(1-t^2)^4\,dt=\|f\|_1.$$ $\endgroup$ May 6 at 1:58
  • $\begingroup$ @ThomasAndrews Thanks for commenting. I am not 100% sure since I found it through Wolfram-Alpha, but it gives the same result under different "forms" as I explained on the question. Even so, following Wolfram-Alpha, it is happening somehow that $\lim_{w \to 0} \hat{f} = \frac{256}{315} \cong \|f\|_1$ (the link is in the comment of Q3)... Is ever weirder since if I split the function on each individual fraction all are singular by themselves, but added they are accurately forming a function that resembles the Gaussian spectrum, except near zero. $\endgroup$
    – Joako
    May 6 at 2:32
  • 2
    $\begingroup$ If evaluated in the straightforward way using computer floating-point arithmetic, the numerator will involve the subtraction of two numbers very close together, which is the very worst situation for precision of floating-point estimates. $\endgroup$
    – aschepler
    May 6 at 2:37
  • 1
    $\begingroup$ this shows the local behavior at 0: wolframalpha.com/… $\endgroup$ May 6 at 3:07
  • 1
    $\begingroup$ @Joako there is a simple reason for what you saw in your edit: for any sufficiently nice function $f$, $\frac1{2\pi}\int_{\mathbb R}\hat f(w)dw = \mathcal F^{-1} \hat f(0) = f(0)$. In your case, $f(0)=g(0)=1$. $\endgroup$ May 10 at 2:34

2 Answers 2

7
$\begingroup$

Computation by hand of the Fourier transform

First, recall that (with your convention for the Fourier transform), $$ \mathcal F[\theta(1-t^2)](w) = \mathcal F \unicode{x1D7D9}_{[-1,1]}(w)= 2 \DeclareMathOperator{\sinc}{sinc}\sinc w. $$

Second, recall that for any polynomial $p(t)$, and any reasonable function (or even distribution) $g$, $$ \mathcal F ( p(t) g)(w) =p(D) \mathcal Fg(w), \qquad D:= \frac1i\frac{d}{dw}.$$ Since $(1-t^2)^4 = 1 - 4t^2 + 6t^4 - 4t^6 +t^8,$ $$ (1-D^2)^4 = 1 +4\frac{d^2}{dw^2} +6\frac{d^4}{dw^4} + 4\frac{d^6}{dw^6} + \frac{d^8}{dw^8},$$ and hence for your function $f(t) = (1-t^2)^4 \unicode{x1D7D9}_{[-1,1]}(t) $, $$\mathcal Ff(w)= 2\sinc w + 8\sinc^{(2)}w + 12 \sinc^{(4)}w+8\sinc^{(6)}w+2\sinc^{(8)}w.$$ Using Wolfram to do the tedious calculations (screenshot) gives the same formula as you, so the answer is right. As $\sinc$ is entire, all its derivatives are, and hence $\mathcal F f$ is as well.

Verifying the behaviour near $0$

Note that in the formula $$ \mathcal Ff(w)=\frac{768 \Big(5w\,(2 w^2-21)\cos(w)+(w^4-45w^2+105) \sin(w)\Big)}{w^9},$$ if we want to see that the numerator is $O(w^9)$, we should (unless we are super lucky) use at least a 8th order Taylor expansion around 0 for the cosine on top, and a 9th order expansion for the sine. This is because of the low order terms in the polynomial coefficients ($-21$ and $+105$, respectively.)

Denote the $n$th partial Taylor expansion around $x_0$ of $f$ at $x$ by $$ T^n f(x_0;x) = \sum_{k=0}^n \frac{f^{(k)}(x_0)}{k!}(x-x_0)^n.$$ Since all derivatives of sine and cosine are uniformly bounded by 1, Taylor's theorem gives, $$ |\cos(w) - T^8 \cos(0;w)| \le \frac{w^9}{9!}, $$ and similarly $$ |\sin(w) - T^9 \sin(0;w)| \le \frac{w^{10}}{10!} ,$$ If we just replace $\cos(w)$ with $T^8 \cos(0;w)$ and $\sin(w)$ with $T^9 \sin(0;w)$ to get $$\tilde f(w):=\frac2{945} (w^4 - 27 w^2 + 384),$$ we incur an error of size at most $$ \frac{768 \Big(5|w|\,(2 w^2\color{red}+21)|w|^9/9! +(w^4\color{red}+45w^2+105) w^{10}/10!\Big)}{|w|^9} = \frac{|w| (1155 + 145 w^2 + w^4)}{4725}. $$ For say $|w|<0.2$, this error is at most $$\frac{1160.8016|w|}{4725}<\frac {|w|}4<\frac1{20}.$$ Therefore, $\mathcal Ff$ stays in the shaded blue region below. The approximation $\tilde f$ is plotted in red. (Desmos link) enter image description here Of course, using more terms will give a more impressive bound, but this is enough to rule out the spiky behavior you saw.

In fact, in the current situation there is extra cancellation, since $T^8 \cos(0;w)=T^9 \cos(0;w)$ and $T^9 \sin(0;w)=T^{10} \sin(0;w)$, giving us a much better error estimate. Redoing the above calculations shows that the error is no more than $(w^2 (w^4 + 155 w^2 + 1260))/51975$ which is $\le w^2/30$ for $w\le 1$. This explains why the approximation is a much better match than the above error plot suggests. The true function $\mathcal Fg$ is plotted in black (for $|x|>0.2$) below: (desmos link)

enter image description here

You have to zoom in quite far in order to even notice that they don't exactly match at $x=0.2$: enter image description here

$\endgroup$
10
  • 1
    $\begingroup$ I don’t see any benefit for using them, but it would not be hard to hard code the correct behaviour near zero in general practice by replacing with a good Taylor series near the origin. If you change the polynomial you get a different combination of derivatives and desmos already struggles with the sixth derivative of sinc. I haven’t watched the video but a Gaussian decays much much faster than 1/w and is practically compactly supported in many circumstances $\endgroup$ May 9 at 2:06
  • 1
    $\begingroup$ I've had to do some analysis around $0$ like this for stuff in my work (different integral transforms, same idea though) because of numeric instabilities due to "divide by zero" issues. I approximate by a Taylor series at $0$ as well. Great answer. $\endgroup$ May 9 at 17:57
  • 1
    $\begingroup$ Sure! It's a mixture of integral transform theory and functional analysis in the vein of quantum mechanics. $\endgroup$ May 10 at 5:55
  • 1
    $\begingroup$ Hi @Joako, I tried to click on the video link but it seems to link to Wikipedia? $\endgroup$ May 13 at 13:38
  • 1
    $\begingroup$ @Joako I have seen your newer question(s) but am stretched for time. I do occasionally think about them and if I have something to say I will comment there. Of course you can always leave me messages in the previously linked chatroom $\endgroup$ Jun 2 at 7:18
3
$\begingroup$

Your $\hat f$ is correct. The graphing programs are not able to give accurate graphs of that $\hat f$ function near $x=0$, at least not from that formula.

The most common form of representing real numbers in computer languages is floating-point arithmetic, where a possible value is an integer times a power of two. The double type lets the integer part use up to 52 bits, so relative errors from operations like multiplication tend to be a few times $2^{-52}$ in magnitude relative to the value. But when addition or subtraction should result in a number much closer to zero, many of the higher-order bits cancel out, and a small proportion of the value bits could be accurate.

In the formula

$$ \hat{f}(w) = \frac{768 \Big(5w\,(2 w^2-21)\cos(w)+(w^4-45w^2+105) \sin(w)\Big)}{w^9} $$

when $w$ is near zero, the addition in the numerator adds a number approximately equal to $-105w$ and a number approximately equal to $105w$, and we expect to get a number approximately equal to $w^9/945$, much much closer to zero. So since plotting programs probably use these double values to be able to get many function values in a reasonable time, they'll be inaccurate when trying to get values near zero.

We can see similar errors asking Wolfram to plot

$$ g(x) = \left(\sin x - x + \frac{x^3}{6} - \frac{x^5}{120}\right)\frac{7!}{x^7} $$

If I convert your $\hat f$ to its Taylor expansion

$$ \hat{f}(w) = 768 \sum_{k=0}^\infty \left( \frac{1}{(2k+1)!} - \frac{10}{(2k+2)!} + \frac{45}{(2k+3)!} - \frac{105}{(2k+4)!} + \frac{105}{(2k+5)!} \right) (-1)^k w^{2k-4} $$

Wolfram still gets the same bumpy graph, possibly by doing symbolic analysis to get back to the sines and cosines. (Take off the word "plot" from the input, and it shows the equivalent formula in terms of sine and cosine, then the same plots again. So at least this tells me I got the series right.)

But then if I limit that Taylor sum to a finite number of terms like $k \leq 20$, instead of asking for the infinite sum, Wolfram can give a more accurate plot, which is indeed smooth near zero. With terms up to $w^{36}$, that should be plenty of accuracy enough as long as $|w|$ isn't huge, and of course a Taylor expansion is at its most accurate near zero.

$\endgroup$
5
  • $\begingroup$ Thanks you very much for the complete answer. When you are saying "near zero" you are meaning "the point where cancellatiom should happen", right?... I plot a displaced version $w \to w-40$ and the rippling singularity is now plot at $w=40$ here... so is not an issue of ploting transforms at the origin. This have never happen to me before. Thanks a lot! $\endgroup$
    – Joako
    May 6 at 6:23
  • 1
    $\begingroup$ Right, this sort of large floating-point inaccuracy happens when the result should be close to zero. $\endgroup$
    – aschepler
    May 6 at 14:16
  • 2
    $\begingroup$ @Joako : The problem is not so much cancellation but the design of the numerical function for sine and cosine. Roughly they are min-max approximations over some interval like $[0,\frac\pi4]$ of $\sin(x)-x$ with odd powers starting from the third and $\cos(x)-1$ with even powers starting from the second. The degrees of the polynomials are chosen so that the error is uniformly bounded by the machine precision over these intervals. This in turn means that the higher-degree coefficients are not exactly the Taylor coefficients, which produces unavoidable error terms in the formula. $\endgroup$ May 10 at 9:11
  • $\begingroup$ @LutzLehmann with "unavoidable" I should interpret that there is no way to find $\|\hat{f}\|_1$ without using approximations? (mysteriously, at least for me, Wolfram-Alpha do find $\int_{\mathbb{R}} \hat{f}\,dw$ and $\|\hat{f}\|_2^2$, and since $\|\hat{g}\|_1 < \infty$ I was expecting to find it but I don't know if it is diverging, or it is because this numerical issue Wolfram-Alfa cannot find it) $\endgroup$
    – Joako
    May 11 at 19:50
  • 1
    $\begingroup$ It is unavoidable if you use the numerical implementations of these functions. A CAS usually works with symbolic representations of the functions. $\endgroup$ May 11 at 20:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.