Unwanted sine modes appearing in convolution using FFT

[pw0908]

Hi All, I don't know if this is an issue with FFTW.jl specifically or something that I've done, but any guidance that could be provide will be appreciated. I've also posted this issue in FastTransforms.jl in case the issue originates from their package: JuliaApproximation/FastTransforms.jl#257

To give a little bit of background, I'm trying to perform the following convolution (nV, r and r' are vectors):

$$n_\mathbf{V}(r) = \int d\mathbf{r}' \rho (r') \delta(|r-r'|-R) \frac{r-r'}{|r-r'|}$$

Where rho(r) is a profile which looks like (I've made sure it's periodic):

I have obtained an analytical expression for the fourier transform of the last two terms in the integral as:

$$\hat{w} = 4\pi i \frac{k}{|k|^3} (\sin(|k|R)-|k|R\cos(|k|R))$$

If I use just regular Float64, the resulting nV profile looks fine:

Except when I zoom in:

While these values are small, in a later part of my code, I need to divide this profile by a small number which makes this 'noise' blow up and affects calculations downstream.

I tried using FastTransforms.jl since it would let me go to higher precision. However, even when I use BigFloat, although it removes the noise, some unphysical oscillations remain:

Im no longer certain that these remaining oscillations are due to floating point errors. I also call these oscillations unphysical as, aside from the nature of this convolution integral, if I perform this convolution integral manually, the oscillations disappear completely (and I'm able to obtain the integral to higher precision).

Am I doing something wrong in the way I'm performing the convolution integral? Thanks!

I've attached an MWE below:

using FFTW, FastTransforms

# Generating the density profile
tanh_prof(x,start,stop,shift,coef) = 1/2*(start-stop)*tanh((x-shift)*coef)+1/2*(start+stop)

ub = 10
lb = -10
z = LinRange(lb,ub,4000)
rho1 = 1e3
rho2 = 1e-12
rho = @. tanh_prof(z,rho1,rho2,(ub/4+3*lb/4),20)*(z<=0) +
         tanh_prof(z,rho2,rho1,(3*ub/4+lb/4),20)*(z>0)

# Obtaining the fourier transform of w_hat
freq = fftfreq(4000,4000/20)
R = 1/2*2pi

weight_hat = @. 0.0 - 4π*im*freq/abs(freq)^3*(sin(abs(freq)*R)-R*abs(freq)*cos(abs(freq)*R)) *(freq != 0.0)

# Performing the convolution integral
rho_hat = fft(rho)
I_hat = @. rho_hat*weight_hat
nV = real(ifft(I_hat))

[stevengj]

You aren't computing a Fourier transform. You're computing a discrete Fourier transform, which has (a) discretization and (b) periodic boundaries. Invariably, one of these two facts explains the differences that people see compared to analytical Fourier transforms and integral convolutions.

[pw0908]

Thanks for your reply @stevengj. Just to clarify a few things:

1. My starting profile is already periodic (or at least, this profile is meant to represent one period of the wave)
2. I'm not sure if this is what you mean by discretization but one other observation I made is that increasing the number of points, even by an order of magnitude, doesn't remove these oscillations. In fact, their amplitude stays the same.

I guess my only question is: is this issue just inherent to doing a discrete fourier transform and, as such, can never be eliminated?

[stevengj]

What I can tell you is that the FFT calculation is correct; it has been tested for decades now, and nothing in your issue leads me to believe that there's anything new here. If you have a bug in your math it is elsewhere, but it's not something I have time to help you with.