DFT derivative property?

Question

Does it have one? The continuous variant does, $f'(t) \rightarrow j \omega F(\omega)$ - but $jkX[k]$ definitely isn't it for DFT.

To find it there must be a useful simplification of $\text{DFT}(x[n] - x[n-1])$, where $x[n] = \sum_{k=0}^{N-1} X[k] e^{j2\pi kn/N}$, while also accounting for the one dropped sample (x[0]) - and I cannot see such a simplification.

---

Note: $x[n] = ...$ above is for the definition of IDFT I use, rather than where the signal derives from. The signal should be assumed the most general possible; not periodic, can be complex or real, but must be finite.

Note2: see my answer for some important caveats to (and direct contradicting of) the accepted answer, and optionally the comments discussions below answers. Thanks to @CedronDawg for the detailed discussion.

Note3: "derivative" = finite difference (for sake of this question).

Answer 1

To get a detailed answer along the lines of what you propose, we need to be careful about the normalization used in discrete Fourier transform (DFT) and inverse discrete Fourier transform (IDFT):

$$\text{DFT: }\quad X[k] = \sum_{n=0}^{N-1} x[n] e^{-j 2\pi kn/N}\tag{1}$$ $$\text{IDFT: }\quad x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j 2 \pi k n / N}\tag{2}$$

Those have a normalization that is directly compatible with fft and ifft from MATLAB, Octave, NumPy and SciPy. The indexes $k$ and $n$ run from $0$ to $N-1$. (MATLAB and Octave have a different indexing convention, $1$ to $N$.) Then:

$$\begin{array}{c}y[n] = x[n] - x[\operatorname{mod}(n-1, N)]\\ \begin{align}\\ Y &= \operatorname{DFT}\big(y\big)\\ &= \operatorname{DFT}(x*[1, -1, 0, 0, \ldots])\\ &= \operatorname{DFT}(x)\times\operatorname{DFT}([1, -1, 0, 0, \ldots])\\ &= X\times\operatorname{DFT}([1, -1, 0, 0, \ldots])\end{align}\end{array}\tag{3}$$ $$\Rightarrow Y[k]= X[k]\left(1 - e^{-j2\pi k/N}\right),\tag{4}$$

where $\operatorname{mod}$ gives the unsigned remainder, for example $\operatorname{mod}(-1, N) = N - 1$, the symbol $*$ denotes length-$N$ circular convolution and $\times$ denotes multiplication, and all sequences are of length $N$.

Circular convolution in the discrete time domain is exactly equivalent to multiplication in the discrete frequency domain, when DFT and IDFT are used to transform the sequences between the domains. See for example Circular Convolution - MIT OpenCourseWare. For "linear" convolution of discrete sequences, there is no such elegant pair of equivalent operations, which makes me think no expression proposed as an answer to the question will ever be fully satisfactory.

Considering input signals that extend to the left of the output range $0 \le n < N$ of IDFT, calculating the backward difference using modulo indexing is conditionally equal to calculating it without it:

$$x[n] - x[\operatorname{mod}(n-1, N)] = x[n] - x[n-1]\quad\text{conditionally}\tag{5}$$

under the condition that you only calculate it for some of the indexes:

$$0 < n < N,\tag{6}$$

or for $0 \le n < N$ under the condition that there's a hint of periodicity in the signal:

$$x[-1] = x[N-1].\tag{7}$$

A sufficient but not necessary condition that satisfies Eq. 7 is that $x$ is $N$-periodic, which is however prohibited by the signal defined as finite in the question. An example of another condition that satisfies Eq. 7 is $x[-1] = x[N-1] = 0$.

@CedronDawg's first answer provides a formula which corrects in the frequency domain the error in Eq. 5 if neither condition is satisfied: $Y[k] = X[k]\left( 1 - e^{-j2\pi k/N} \right) - x[-1] + x[N-1]$. As an alternative derivation, in length-$N$ time domain the correction is an impulse:

$$\begin{align}&\big[x[0] - x[−1] - \big(x[0] - x[N-1]\big),\, 0,\, 0,\, \ldots\big]\\ = &\big[x[N-1] - x[-1],\, 0,\, 0,\, \ldots\big],\end{align}\tag{8}$$

which in frequency domain is a constant (see DFT Pairs and Properties: pair row 2, property linearity):

$$x[N-1] - x[-1],\tag{9}$$

to be added to all elements of $Y$ calculated by Eq. 4.

For a general $x$, the condition of Eq. 6 for Eq. 5 enables to use a length $N+1$ DFT and IFT to calculate the backward difference, by shifting the input to the DFT one step to the right, and finally by shifting the output from IDFT one step to the left. With forward difference $x[n + 1] - x[n]$, the shift would not be necessary, and I think this matches your discarding of the 0th sample. For a circular convolution implementation of convolution by a finite sequence, using a longer transform is a common trick to avoid the circular effects in a sufficiently large part of the output of the IDFT. For then obtaining the DFT of a partial IDFT output, in particular a one shorter, I don't think there is any shortcut.

We could perhaps re-express the question as: What is the DFT of the length $N-1$ forward difference of $x$ of length $N$, given $x$ and $X_{N-1} = \operatorname{DFT}(x_{N-1})$, a length $N-1$ DFT of the partial sequence $x_{N-1} = \big[x[0], x[1], \ldots, x[N-2]\big]$? Analogously to Eq. 4 we have:

$$y_{N-1} = x_{N-1}*[-1, 0, 0, \ldots, 0, 0, 1]\tag{10.1}$$ $$\Leftrightarrow Y_{N-1} = X_{N-1}\times\operatorname{DFT}([-1, 0, 0, \ldots, 0, 0, 1])\tag{10.2}$$ $$\Rightarrow Y_{N-1}[k] = X_{N-1}[k]\big(e^{j2\pi k / (N - 1)} - 1\big),\tag{10.3}$$

where each sequence is of length $N-1$. The desired forward difference $f_{N-1}$ is:

$$f_{N-1} = \big[x[1] - x[0],\, x[2] - x[1],\, \ldots,\, x[N-1] - x[N-2]\big].\tag{11}$$

Eq. 10.1 can be expanded to:

$$y_{N-1} = \big[x[1] - x[0],\, x[2] - x[1],\, \ldots,\, x[0] - x[N-2]\big].\tag{12}$$

By comparing Eqs. 11 and 12, it can be seen that:

$$f_{N-1} = y_{N-1} + \big[\ldots,\, 0,\, 0,\, x[N-1] - x[0]\big],\tag{13}$$

where the sequence in brackets is of length $N-1$. Taking the DFT of both sides of Eq. 13 and applying Eq. 10.3 gives the answer:

$$\begin{align}F_{N-1}[k] &= Y_{N-1}[k] + (x[N-1] - x[0])e^{j2\pi k/(N-1)}\\ &=X_{N-1}[k]\big(e^{j2\pi k / (N - 1)} - 1\big) + (x[N-1] - x[0])e^{j2\pi k/(N-1)}.\end{align}\tag{14}$$

This is the length $N-1$ DFT of the length $N-1$ forward difference of $x$ of length $N$.

Alternatively, you might be interested in the derivative of the band-limited signal represented by the samples.

Answer 2

[Edit: This answer is based on x[n] being the input signal (likely sampled from a continuous differentiable function), which was my misinterpretion of the OP's question, and that what was being sought was a transform from the spectrum of the input signal to the spectrum of the linear first difference. If one desires the spectrum of the circular first difference within the analysis frame there is no need for the slope term. Furthermore, it can clearly be seen from the formula that the DC bin of the circluar difference spectrum will be zero.]

This is a different approach than Olli's doing it straight from the definition and not assuming a N periodic signal.

$$ \begin{aligned} Y[k] &= \text{DFT}(x[n] - x[n-1]) \\ &= \text{DFT}(x[n]) - \text{DFT}(x[n-1]) \\ &= X[k] - \text{DFT}(x[n-1]) \\ &= X[k] - \sum_{n=0}^{N-1} x[n-1] e^{-i\frac{2\pi}{N}kn } \\ &= X[k] - \sum_{m=-1}^{N-2} x[m] e^{-i\frac{2\pi}{N}k(m+1)} \\ &= X[k] - e^{-i\frac{2\pi}{N}k} \sum_{m=-1}^{N-2} x[m] e^{-i\frac{2\pi}{N}km} \\ &= X[k] - e^{-i\frac{2\pi}{N}k}\left[ \sum_{m=0}^{N-1} x[m] e^{-i\frac{2\pi}{N}km} + x[-1]e^{i\frac{2\pi}{N}k} - x[N-1] e^{-i\frac{2\pi}{N}k(N-1)} \right] \\ &= X[k] - e^{-i\frac{2\pi}{N}k}\left[ X[k] + x[-1] e^{i\frac{2\pi}{N}k} - x[N-1] e^{-i\frac{2\pi}{N}k(N-1)} \right] \\ &= X[k]\left( 1 - e^{-i\frac{2\pi}{N}k} \right) - \left[ x[-1] - x[N-1] e^{-i\frac{2\pi}{N}k(N-1)}e^{-i\frac{2\pi}{N}k} \right] \\ &= X[k]\left( 1 - e^{-i\frac{2\pi}{N}k} \right) - x[-1] + x[N-1] \\ \end{aligned} $$

For a N periodic signal $x[-1] = x[N-1] $ so the result becomes:

$$ Y[k] = X[k]\left( 1 - e^{-i\frac{2\pi}{N}k} \right) $$

Considering the normalization is important. I used the conventional unnormalized forward DFT definition. Using the (more proper,IMO) 1/N normalization, as implied by the OP's definition of the inverse DFT, the expression becomes:

$$ Y[k] = X[k]\left( 1 - e^{-i\frac{2\pi}{N}k} \right) + \frac{x[N-1]-x[-1]}{N} $$

The last term now clearly becomes a slope calculation.

---

In response to OverLordGoldDragon's comment. Using the conventional normalization:

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{i \frac{2\pi}{N} nk } $$

Let's make it continuous.

$$ x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{i \frac{2\pi}{N} nk } $$

Take the derivative in respect to n.

$$ \begin{aligned} \frac{dx}{dn} (n) &= \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{i \frac{2\pi}{N} nk } \left(i \frac{2\pi}{N} k \right) \\ &= \frac{2\pi}{N^2} \sum_{k=0}^{N-1} i k X[k] e^{i \frac{2\pi}{N} nk } \\ \end{aligned} $$

So, your initial assertion isn't quite correct, it does apply to the DFT for the continuous derivative.

By the chain rule:

$$ \frac{dx}{dt} = \frac{dx}{dn} \cdot \frac{dn}{dt} $$

The differential you are using for an approximation for the derivative is not the best one as it has a half sample shift included. Generally you would prefer $(x[n+1]-x[n-1])/2$. This stays centered and doesn't "amplify noise" in the higher frequency range.

The comment is a little off too. $$ \begin{aligned} \left( 1 - e^{-i\frac{2\pi}{N}k} \right) &= e^{-i\frac{\pi}{N}k} \left( e^{i\frac{\pi}{N}k} - e^{-i\frac{\pi}{N}k} \right) \\ &= e^{-i\frac{\pi}{N} k} 2i \sin\left( \frac{\pi}{N}k \right) \\ &= e^{-i\frac{\pi}{N} k} 2 e^{i \frac{\pi}{2} } \sin\left( \frac{\pi}{N}k \right) \\ &= e^{-i\left( \frac{\pi}{N} k - \frac{\pi}{2} \right) } 2 \sin\left( \frac{\pi}{N}k \right) \\ \end{aligned} $$

so

$$ \left| \left( 1 - e^{-i\frac{2\pi}{N}k} \right) \right| = 2 \sin\left( \frac{\pi}{N}k \right) $$

and

$$ \arg \left( 1 - e^{-i\frac{2\pi}{N}k} \right) = -\frac{\pi}{N} k + \frac{\pi}{2} = - \frac{\pi}{2} \left( 1 - \frac{2k}{N} \right) $$

Switching gears, note that for $k=0$

$$ Y[0] = \frac{x[N-1]-x[-1]}{N} $$

So the DC component of your differential is the average slope over your frame, as expected.

I don't read text books, so I can't really address that. Interesting is not always synonymous with useful.

I'm self taught, then I do online research to confirm my findings. What is a wonder to me is that my frequency formulas, which are the first to achieve exactness, aren't in the curriculum yet either.

---

Puzzle solved. Duh.

Unless you center around zero, so you are using $k=-1$ instead of $k=N-1$, the derivative will be of the higher frequency interpolation (equivalent to the derivative of the DTFT at that point).

This is similar to the "Fluffy Cloud" case here: How to get Fourier coefficients to draw any shape using DFT?

---

I think part of my confusion is that you are using (upon rereading) the extended inverse DFT as x[n], thus x[-1] can be calculated and will match x[N-1], while in your code you generate a signal x[n] from scratch.

Define

$$ y[n] = x[n] - x[n-1] $$

and

$$ Y[k] = DFT(y[n]) $$

Your code "drops a sample", while Olli's and my answers employ x[-1]. Olli's answer assumes periodicity (go ahead accept it), mine doesn't (thinking you were working with a raw signal). I don't think it is appropriate to "drop a sample" as the np.diff call does. If you do, you should reframe the problem on the domain of 1 to N-1 as being 0 to M-1, where M = N-1, then you have the equivalent situation as not having dropped a sample.

When your "drop a sample" it changes the DFT definitions invalidating my first two lines.

Be a bit patient, and I'll clarify the puzzle remark.

---

Without a lot of explanation, here is the "puzzle solved". Note, that if x were to be interpolated using the upper k values as positive frequencies, there would be a whole lot of oscillations between the plotted points. Thus, so would Y.

Rescaling is muddled by the $ 2\pi $ factor on t and the $\frac{dn}{dt}$ factor, so I didn't bother cluttering the code as it isn't salient to the issue at hand.

import numpy as np
import matplotlib.pyplot as plt

#=========================================================
def main():

        N = 128
        
        t = np.linspace( 0, 1, N, False )
        x = np.cos( 2 * np.pi * t )
        
        X = np.fft.fft( x )

        plt.plot( x )
        plt.show()
        
        Y = np.zeros( N, dtype='complex' )
        
        for k in range( N ):
          Y[k] = X[k] * 1j * k

        y = np.fft.ifft( Y )
        
        plt.plot( y.real )
        plt.plot( y.imag )
        plt.show()
        

        Z = np.zeros( N, dtype='complex' )
        
        H = N >> 1
        
        for k in range( -H, H ):
          if k >= 0:
             Z[k]   = X[k]   * 1j * k
          else:
             Z[k+N] = X[k+N] * 1j * k

        z = np.fft.ifft( Z )
        
        plt.plot( z.real )
        plt.plot( z.imag )
        plt.show()


#=========================================================
main()

Here is the last plot:

Here is a test program for the formula:

import numpy as np

#=============================================================================
def main():

        L  = 10000              # Length of Signals
        F  = 100                # Frame location 

        N  = 16                 # Frame Size = DFT Sample Count
        
        P = np.zeros( L )       # Position
        D = np.zeros( L )       # Difference
        
        for n in range( L ):
          P[n] = 1.2 + 0.3 * n + 0.045 * n * n

        for n in range( 1, L ):
          D[n] = P[n] - P[n-1]
          
        x = P[F:F+N]  
        y = D[F:F+N]
        
        X = np.fft.fft( x ) / N
        Y = np.fft.fft( y ) / N
        
        Z = np.zeros( N, dtype=complex )
        
        C = ( x[N-1] - P[F-1] ) / N

        for k in range( N ):
          Z[k] = X[k] * ( 1 - np.exp( -1j * ( 2.0 * np.pi / N ) * k ) ) + C
        
        for n in range( N ):
          print( "%2d %10.6f %10.6f   %10.6f %10.6f" % \
               ( n, Y[n].real, Y[n].imag, Z[n].real, Z[n].imag ) )
        

#=============================================================================
main()

Here are the results:

 0   9.930000   0.000000     9.930000   0.000000
 1  -0.045000   0.226230    -0.045000   0.226230
 2  -0.045000   0.108640    -0.045000   0.108640
 3  -0.045000   0.067347    -0.045000   0.067347
 4  -0.045000   0.045000    -0.045000   0.045000
 5  -0.045000   0.030068    -0.045000   0.030068
 6  -0.045000   0.018640    -0.045000   0.018640
 7  -0.045000   0.008951    -0.045000   0.008951
 8  -0.045000   0.000000    -0.045000  -0.000000
 9  -0.045000  -0.008951    -0.045000  -0.008951
10  -0.045000  -0.018640    -0.045000  -0.018640
11  -0.045000  -0.030068    -0.045000  -0.030068
12  -0.045000  -0.045000    -0.045000  -0.045000
13  -0.045000  -0.067347    -0.045000  -0.067347
14  -0.045000  -0.108640    -0.045000  -0.108640
15  -0.045000  -0.226230    -0.045000  -0.226230

Answer 3

Response to RB-J's comments under my other answer:

"what does a DFT periodicity denier do with x[−1]?"

"what if you multiply the DFT of x[n] (which we call "X[k]") with ej2πdk/N (where d is an integer)? you will find that the DFT very much cares about periodicity. – "

And similar from Making the units of a analytically calculated PSD consistent with the units of an FFT

"//the DFT does not consider what the signal is outside the frame// ... never? ---okay @CedronDawg, what does the DFT "consider" when you multiply X[k] with ej2πdk/N where d∈Z? (all of the X[k] for 0≤k<N.) – robert bristow-johnson 4 hours ago "

Okay, let's start with a signal with lots and lots of samples. Suppose I want to analyze a section of it with N samples in it. To do so, I will establish an analysis frame. This means re-indexing the samples so that the first sample of my analysis frame is at index 0. Since I intend to take a DFT on this section, I can also call it a DFT frame.

This is not the same thing as "windowing". I could also define a rectangle window function having ones for all the samples in my section, and zeroes everywhere else. The window function is indexed on the underlying indexing, it does not establish a new indexing frame.

If I wanted to, I could also define a window function relative to my analysis frame. A rectangle window spanning my frame would be inert/moot in terms of taking the DFT. You can consider the DFT definition as having an implicit window like this, but it is not part of the definition. If it were the unnormalized definition would be:

$$ X[k] = \sum_{n=0}^{N-1} 1 \cdot x[n] \cdot e^{-i \frac{2\pi}{N} kn } $$

So, if the section is not at the start of the signal, $x[-1]$ would refer to the sample just prior to my analysis frame. If the section is at the start of the signal then the correct answer would be that $x[-1]$ is undefined. If one had to make an assumption about it, the most logical answer would be a value of zero.

Now, suppose we take the DFT of this section which yields $ X[k] $ values for $ 0 <= k < N $. Values outside this range can also be evaluated by the defintion and will result in a periodic spectrum with a period of N. This may not be the fundamental period, but that is irrelevant to this discussion.

We are now at the starting point of the OP's question (something I misunderstood at first).

Given the values of $X[k]$ we can apply the inverse DFT and reproduce the $x[n]$ values for $ 0 <= n < N $, the samples within the analysis frame. As with the spectrum, the formula also can be evaluated for $n$ values outside the analysis frame. Likewise, the resulting values will form a periodic sequence with period N. There is absolutely no reason to assume, assert, or define that these extended values should match the original signal values outside the analysis frame. If you wish to label that as denial, well ....

It is also impossible, given just the $X[k]$ values to determine what the original $x[-1]$ value is. The value obtained by extending the inverse DFT will be the same as $x[N-1]$. Again, there is no foundation to stipulate, assume, assert, maintain, or proclaim that this matches the original signal value.

Suppose now that I multiply/twist up/modulate the spectrum using a factor of $ e^{i \frac{2\pi}{N} dk } $. Okay, let's do the math.

$$ Z[k] = X[k] \cdot e^{i \frac{2\pi}{N} dk } $$

Take the inverse DFT:

$$ \begin{aligned} z[n] &= \frac{1}{N} \sum_{k=0}^{N-1} Z[k] \cdot e^{ i \frac{2\pi}{N} nk } \\ &= \frac{1}{N} \sum_{k=0}^{N-1} X[k] \cdot e^{i \frac{2\pi}{N} dk } \cdot e^{ i \frac{2\pi}{N} nk } \\ &= \frac{1}{N} \sum_{k=0}^{N-1} X[k] \cdot e^{ i \frac{2\pi}{N} (n+d)k } \\ \end{aligned} $$

Let $ m = n + d $ so $ n = m - d $.

$$ z[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] \cdot e^{ i \frac{2\pi}{N} mk } = x[m \text{ mod } N] = x[ ( n + d ) \text{ mod } N] $$

The end result is that the samples in the analysis frame have been rotated by $d$ samples.

As above, $z[n]$ can be extended using the inverse DFT definition outside the analysis frame forming an N periodic sequence. There is even less rationale to assume that those values will match the original signal.

The OP's question was suppose that $ y[n] = x[n]-x[n-1] $, can you find $Y[k]$ from $X[k]$ without going through the taking the inverse DFT, taking the difference, and then taking the DFT.

The OP stated that $x[n]$ came from the inverse DFT of $X[k]$, whereas I mistakenly assumed the original $x[n]$ was known.

Then answer was given elegantly by Olli. And derived from definition in my answer. The condition that $x[-1]= x[N-1]$ is met for the extended reconstructed signal.

Thus spoke this periodicity denier.

Answer 4

Complementing, and based on C. Dawg's answer, discarding the slope addition, the effect on magnitude and phase are

$$ \begin{aligned} |X[k]| & \rightarrow M|X[k]| \\ \angle{(X[k])} & \rightarrow \phi + \angle{(X[k])} \\ \end{aligned} $$

where

$$ \begin{aligned} M &= 2 \left| \sin{\left( \pi \frac{k}{N} \right)} \right| \\ \phi &= \frac{\pi}{N}(k\ \text{mod}\ N- N/2) \cdot \lceil k\ \text{mod}\ N \rceil \end{aligned} $$

The $\lceil k\ \text{mod}\ N \rceil$ sets $\phi$ to 0 when $k$ is a multiple of $N$, using the convention $\angle(0 + 0i)=0$. Graphically,

Note that the unwrapped $\phi$ is a straight line, so the time-domain effect is a time shift.

---

iDFT: indices should run from 1 to N - 1, as the finite difference drops a sample. Also, the slope addition is redundant; the inverse DFT extends the original signal N-periodicically, so it's zeroed.

Python implementation below; tested with random normal noise, which is a sort of "un-nicest" signal - mean absolute error is 1e-16, which is simply float error.

def d_idft(coef):
    N = len(coef)
    coef = coef * (1 - np.exp(-1j * 2 * np.pi / N * np.arange(N)))
    return np.fft.ifft(coef)[1:]

---

APPENDUM: below are my comments, trimmed from discussions below answers, summarizing important points/caveats, and contradicting the accepted answer. Didn't include other speakers as it'd get really long, but shown excerpts should hint at what's being responded to.

---

Slope term:

"discrete derivative meaningless" - this isn't about the discrete derivative, but about x[n] - x[n-1], whatever interpretation it may hold. Sometimes it's as good as a derivative, other times it's exact in the sense of undoing cumsum, yet other times virtually useless, but point is it's some time-domain transformation whose frequency-domain equivalent we seek

As to your discrete derivation; since there isn't an x[-1] to begin with, the only alternative is the inverse, which equals x[N-1], so again it drops.

It's an overcomplete representation. To encode an N-1 point derivative, it requires N+1 points of data (the original N coefficients plus x[-1]). The DFT is a complete encoding; we can invert it and differentiate in time domain without ever needing x[-1], so same must be possible working purely in the frequency domain.

Consider a counterexample; suppose we don't take x[-1] to come from inversion; then it must be of the original signal, which we framed, but where x[-1] exists. Suppose x[0] to x[N-1] are all zeros, and x[-1] = 100^100 ...

---

"DFT assumes the input is periodic"

I said the inversion, when extended, is periodic; the DTF describes only the portion of the original signal which was fed - not more, not less. We can extend in time domain to analyze aliasing etc., but can't make any description about the original signal outside of the frame.

The root of the debate is a lack of sufficient definitions, i.e. what "assumes" even means; you and I see it as non-fundamental to the transform itself, yet we also agree some contexts / operations demand periodicity.

---

Continuous derivation ($x(n)$)

An inverse "inverts" the forward transformation, and a forward transform of a continuous signal would have $N=\infty$. DTFT: the encoding inserts as much information (i.e. ω) as is decoded (inverted), and can thus be discretized arbitrarily. To contrary, the DFT of a signal demands as many coefficients as the number of points in the signal - whereas your inverse assumes a finite (or simply different) N, hence not decoding what was encoded.

It absurdly implies ability for infinite compression, as you're restoring an infinite x(n) from N finite coefficients.

But what if we simply discretize, $x'(n)\rightarrow x'[n]$? Then the compression argument falls. Your expression handled a linear chirp fairly well in terms of shape, comparing with np.diff, and since latter indeed isn't a "true" derivative, it's not obvious which is "more correct". There's math work to be done to justify continuous-differentiating and then discretizing again, but maybe it can be done after all.

Edit 1/10/2020: indeed, it can be done, and is interpreted as a discretized (not discrete) derivative of the trigonometric interpolation of $x[n]$ (not of $x(n)$), and the 'nicer' such interpolation fits $x[n]$, the more accurate the derivative. Stranger yet, discretized FT and DFT can work together sometimes.

However, Cedron's derivation still normalizes wrongly (but his code is correct); strangely, $x'[n]$ is found via discrete iFT (iDFT) of discretized FT derivative, so the extra $2\pi / N$ doesn't belong. ... or maybe it belongs under a different set of assumptions, but none I've yet come across.

Answer 5

WINNER: Olli's solution:

$$ F_{N-1}[k] = X_{N-1}[k]\big(e^{j2\pi k / (N - 1)} - 1\big) + (x[N-1] - x[0])e^{j2\pi k/(N-1)} $$

Code + Demo:

def dft(x):
    return np.fft.fft(x[:-1])

def d_idft(coef, x):
    M = len(x) - 1  # N - 1
    exp = np.exp(1j * (2 * np.pi / M) * np.arange(M))
    coef = coef * (exp - 1) + (x[-1] - x[0]) * exp
    return np.fft.ifft(coef)

Though, problem not entirely solved. All properties 1 & 2 met, negligible reconstruction loss, and very small difference in coefficients computed directly from DFT(DIFF(x)) and this method. "Very small" as in 1e-12, which might be FFT's float error. -- Comparison code.

1. Need a relation in terms of $X_N=\text{FFT}(x)$, not $X_{N−1}=\text{FFT}(x[\text{:-1}])$
2. Cannot require the original signal as input to compute $F_{N−1}$ (irrecoverable from $X_{N−1}$).

Since Olli worked out an expression with 1 fewer coefficient ($X_{N−1}$), in a way a harder problem was solved, so a workaround to (1) and (2) seems possible.

---

This answer aims to further refine the problem, test proposed solutions, and explicate flaws in other answers to guide a resolution.

Slope term: it's needed, whether $x$ is periodic or not. Without it, $Y[0] \neq \sum_{n=1}^{N-1} (x[n] - x[n-1])$. It somewhat ties to the Fundamental Theorem of Calculus; if the finite difference is the derivative, then $x[n]$ is the antiderivative of the function that's the finite difference, so the sum for $Y[0]$ is quite simply $x[N-1]-x[0]$ -- [$\int_a^b f'(x)dx=f(b)-f(a)$].

More precisely, $x[n-1]-x[0]$ equals the net finite difference on the interval; try to force this value to be anything else without changing the endpoints - you won't.

---

Answer properties - what the resultant expression must satisfy:

1. len(Y) == N - 1; the finite difference uses 1 fewer term than x
2. $Y[0] = \sum_{n=1}^{N-1} (x[n] - x[n-1])$
3. $Y$ is expressed in terms of $X$, $j$, $n$, $N$, and maybe its own index - no other parameter; $x$ is fine as long as within $[0, N-1]$, as that's recoverable from $X$.

---

Properties justification:

Olli's answer suggests a shift with a N+1-point DFT&IDFT; this is invalid, and so is any attempt at not dropping a sample. It may prove useful in some applications, but is fundamentally flawed and hence not general-purpose. Here's why.

Recall $+C$ when integrating? Same story: differentiation loses information. Let $x'[n] = x[n] - x[n-1]$. We cannot restore $x[n]$ from $x'[n]$, only its shape; the DC offset is missing. If we know any one original sample, $x[h]$, we can restore the whole signal from $x'[n]$, since we know $x[N-1]-x[0] = \sum x'[n]$, and we'll find $x[n-1]$ via cumulative sum starting at $x[h]$. It's analogous to an Initial Value Problem.

---

Cedron's result,

$$ Y[k] = X[k]\left( 1 - e^{-i\frac{2\pi}{N}k} \right) + (x[N-1]-x[-1]), $$

meets none of the properties, but it clearly steps in the right direction.

---

Cedron's result + mod, tested (mod = $x[0]$ instead of $x[-1]$):

The three substantially different magnitude & phase spectra all invert-transform to almost exactly the same signal. Any error seems attributable to float imprecision. Green here is the "ground truth". -- Code.

What if we modify d_idft per above rationale? i.e., drop coef[0], the DC term:

def d_idft(coef, x, slope=True):
    M = len(coef) - 1
    coef = coef[1:] * (1 - np.exp(-1j * 2 * np.pi / M * np.arange(M)))
    coef += (x[-1] - x[0])
    return np.fft.ifft(coef)

The error grows significant - no good. Also note I pass x as input for convenience only.

---

If what we have already works, why look further? (i.e. MAE within float64) - because it's simply wrong, and just because I haven't found a signal with significant MAE, doesn't mean it doesn't exist, or that there aren't graver implications in practice.

---

Test script: you know your solution works if it passes this script.