How to fix shift and scaling errors in an FFT-powered differentiation program?

Question

So I am trying, like many others before me, to sift through all the celestial alignments of phase shifts and normalization coefficients necessary to make differentiation via Fourier transform work.

I am trying to use as little code as possible, making total use of numpy functionality to keep the code clean while using its array operations, which, I was told, are faster than explicit python loops.

From line 106 of the table of Functional relationships, One-dimensional, in the section about Important Fourier Transforms of the Wikipedia article on Fourier Transforms, I get that the Fourier transform has this property where differentiation in real space is equivalent to multiplication in frequency space such that d^nf/dx^n = ifft[F*(i*w)^n], where F is the fourier transform of f.

Looking for a more solid reference to define the what the w would be in d^nf/dx^n = ifft[F*(i*w)^n], I found a PDF paper on the mathematics of differentiating using DFT. In the paper, the author discusses that there are infinitely many of the so-called trigonometric interpolations for a given function and shows a way to obtain to "less wiggly" solution, which is unique.

So I used his definitions and wrote a small python code to do the differentiations. The code I came up is this one (pay attention to the last two lines, since they will change)

N = 51
a = 0
b = 2*np.pi
X  = np.linspace(a,b, N)

L = abs(b-a)
TwoPiOverL = 2.0*np.pi/L
freqs = np.fft.fftfreq(N,1./N)

# THESE COME FROM [1] #
D1 = TwoPiOverL*freqs*1j     
D2 = -(TwoPiOverL*freqs)**2
#=========================#

S = np.sin(X)
fftS=np.fft.fft(S)
C = np.cos(X)

# COMPUTING THE DERIVATIVES #
Y1 = np.fft.ifft(D1*fftS).real   
Y2 = np.fft.ifft(D2*fftS).real

In the figures of the graphs generated, the notation ~d/dx sin means the derivative of the sine function evaluated according to the Fourier Transform multiplication, instead of the analytical expression.

Test 1

First Test, using only the definitions found in the PDF:

So the first test was a mess. The derivatives are well-interpolated around the middle of the domain, but the end-points are crazy. In the PDF itself, at the beginning of Section 6, the author says:

All of the above applies to spectral differentiation of periodic functions y(x). If the function is non-periodic, artifacts will appear in the derivatives due to the implicit discontinuities at the endpoints.

Sounds about right, and it means that I am screwing up somewhere, because if the entire paper is written discussing periodic signals, and if my signal happen to fall in the non-periodic class, then there is room for improvement if I can use a periodic signal.

But I am using DFT, which means I'm applying a numerical procedure to an array. Discrete data that carries no information about periodicity. So then I thought:

Well, how can I tell my Fourier Differentiator that the signal is Periodic? What does that entail, in code? Where, in the code, must I change in order to convey that information to my program?

What seemed like an answer came from the stackoverflow. In this post about Numerical Differentiation using FFT (yes, redundant, but I'll get there), @dkato gives instructions that I have interpreted as changing the two last lines in my code from

Y1 = np.fft.ifft(D1*fftS).real   
Y2 = np.fft.ifft(D2*fftS).real

to

Y1 = np.fft.ifft(np.exp(D1)*fftS).real   
Y2 = np.fft.ifft(np.exp(D2)*fftS).real

So I did that, and then tested what I would get.

Test 2

Second Test, using @dkato corrections:

As seen from the figure above, by using the corrections provided by @dkato, the end-point wiggling goes away, but the estimation of the derivative is still pretty far off. I am using 51 datapoints so I am guessing that this isn't due to inaccuracies with the expansion. It seems like the way I compute the derivatives is incorrect.

The source of the most loss in accuracy between the analytical first derivative and its Fourier Transform estimation seem to be a shift to the right, whereas the second derivative is, not only shifted but also badly scaled. Funny enough, however, when we use

Y2 = np.fft.ifft(-np.exp(D2)*fftS).real

to evaluate the second derivative we see that all of the error due to shifting goes away, leaving only the scaling. I don't know if that is a property of general FFT-powered second derivatives or if it that is just because I am using a Sine function as my y(x), specifically.

The Point, at last

I am happy that the weird, ugly wiggling at the endpoints disappeared, meaning that I must have, somehow, told my Fourier Differentiator that I am working with periodic signals. But these strange shifting and scalings of my differential estimations are still bothersome, and I can't for the life of my put my finger on what is wrong. As you might have guessed, in part because I have little experience playing with all the moving parts in a DFT.

In this capacity, I have two questions for the community if you would be as kind as to assist me with this problem.

1) If it is true that I have somehow told my fourier differentiator that I am using a periodic function, how exactly did I do it? What is the reasoning behind the code that works and the code that does not?

and

2) What is the source of these weird scale/shift incongruities between the analytical and estimated derivatives? How can I fix them?

Answer 1

Let's first answer the crux of the matter:

What is the source of these weird scale/shift incongruities between the analytical and estimated derivatives? How can I fix them?

The problem starts from the following line:

X  = np.linspace(a,b, N)

By default numpy.linspace includes the endpoint, which introduces a discontinuity in the computed periodic function extension since the endpoints are repeated.

To illustrate how this could be a problem, you may look at the computed function sin(X) for N=4:

sin(0)
sin(2*np.pi/3)
sin(4*np.pi/3)
sin(2*np.pi)

for which the phase increments by 2*np.pi/3 at each subsequent value. Repeating these values to get the periodic extension with period N=4, would yield

sin(0)
sin(2*np.pi/3)
sin(4*np.pi/3)
sin(2*np.pi)
sin(0)         # == sin(2*np.pi)
sin(2*np.pi/3) # == sin(2*np.pi + 2*np.pi/3)
sin(4*np.pi/3) # == sin(2*np.pi + 4*np.pi/3)
sin(2*np.pi)   # == sin(2*np.pi + 2*np.pi)

Where you could see that the phase increments by 2*np.pi/3 at each step except at the block boundary where it stays the same. This small discontinuity makes it much harder to get an good approximation with the trigonometric interpolation basis. The resulting interpolation then contains significant oscillations which are further amplified by the derivation operation.

To get rid of the problem you should exclude the endpoint with (and otherwise leaving your initial solution intact):

X  = np.linspace(a,b, N, endpoint=False)

As far as the solution in @dkato's post, it looks like it only implements a time shift operation. For a sinusoidal signal that may look not too bad, but it has nothing to do with differentiation.

On a final note, to answer the part

If it is true that I have somehow told my fourier differentiator that I am using a periodic function, how exactly did I do it?

Yes it is true, and you did it the moment you started using the Discrete Fourier Transform (and its FFT implementation).