2
$\begingroup$

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).

$\endgroup$
42
  • 1
    $\begingroup$ strictly speaking, there is no derivative operator applied to discrete sequences. there are no differential equations involving discrete sequences, but we do approximate (or emulate) differential equations with difference equations. that said, your "derivative" could be meaningfully replaced with a difference but there are issues involved in doing so. one of the issues is an apparent offset of $\frac12$ sample when the discrete samples are cast back into a continuous function using the reconstruction theorem. $\endgroup$ – robert bristow-johnson Sep 11 '20 at 2:44
  • 1
    $\begingroup$ @OlliNiemitalo well, I was already pretty sure, even without reading, that neither of the answers were trying to prove the unexisting derivative property of DFT, but then the only real problem to be solved in this question is to find the correct question to ask for what is being answered below... :-))) $\endgroup$ – Fat32 Sep 12 '20 at 10:59
  • 1
    $\begingroup$ @OlliNiemitalo Finite difference demands dropping a sample; it operates on two points, starting from x[0], ending with x[N-1], which can only produce N-1 points, unless we impute one somehow. I say drop x[0] since that's what Python's numpy.diff does. $\endgroup$ – OverLordGoldDragon Sep 12 '20 at 11:20
  • 1
    $\begingroup$ @CedronDawg I have to agree with Robert here. You don't need to introduce new fancy math just to show that the DFT implies circular first difference, instead of the linear one. And to get the effect of linear first difference, all you have to do is to append one zero to the end of x[n], and use N+1 length DFT. And btw, DFT = DFS. It's nothing new. And it's circular all the way. But the link between a linear sequence $x[n]$ and its circular (periodic) extension $\tilde{x}[n]$ should be made very clear at the beginning so as to allow readers correctly interpret what you are doing.. $\endgroup$ – Fat32 Sep 12 '20 at 12:53
  • 1
    $\begingroup$ @CedronDawg I said fancy, and not fake :-) so please don't resent back ;-). Things are all about definitions and not real math issues here. In the context of DFT, $x[-1]$ simply refers to $x[N-1]$ due to $x[n]$ being circular (periodic) sequence. If you want to make a reference to a linear operation on the supplied data of $x[n]$ treated as an aperiodic finite length sequence,you can do it in a number of ways,but the most elegant is to use the concept of circular vs linear adjustment on the DFT size N. And finally DFT is DFS and there is no need to re-define it I guess... ;-) $\endgroup$ – Fat32 Sep 12 '20 at 14:40
6
$\begingroup$

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.

$\endgroup$
34
  • $\begingroup$ Seem to have nailed it - I'll run some checks to confirm this works in my context. $\endgroup$ – OverLordGoldDragon Sep 4 '20 at 12:39
  • $\begingroup$ A question as simple as DFT{x[n] - x[n-1]} should be answered as simple as $X[k] (1 - e^{-j\frac{2\pi}{N}k})$, and that arguments are modulo-N.;i.e., $x[n-1] = x[(n-1)_N]$ and that $x[-1] = x[N-1]$. Now, I understand that the OP wants to get a linear first difference instead of circular one, and then the answer should only be modified by adding that $M=N+1$ length cicular difference will be equal to the linear difference. Just as the case of liner vs circular convolution. Nothing else needs to be added. But we all tend to improve already perfect answers bu putting more work of us.;-) $\endgroup$ – Fat32 Sep 12 '20 at 12:42
  • $\begingroup$ Haven't studied circular in detail yet, can't comment there - but the end-result still fails properties 1 & 2: len(Y) == N, and Y[0] != sum(x[n] - x[n - 1]). Prop 2 must hold whether or not $x$ is periodic - but in your answer, Y[0] == 0. $\endgroup$ – OverLordGoldDragon Sep 13 '20 at 11:19
  • 1
    $\begingroup$ Err, a small typo in my code. I have a grand announcement... It WORKS! Excellent work, Olli. $\endgroup$ – OverLordGoldDragon Sep 15 '20 at 21:54
  • 1
    $\begingroup$ I don't think there is a shortcut to "directly relate the frequency spectrum of x with that of DIFF(x)", as you express it. The main computation would amount to taking the IDFT and then a one-shorter DFT, of which I say the same thing in the answer. That computation would be equivalent to band-limited resampling, which is a quite thoroughly researched area of digital signal processing, and I doubt there is anything new to be found there, math or computation-wise, that would also be applicable to the present task. $\endgroup$ – Olli Niemitalo Sep 16 '20 at 5:55
4
$\begingroup$

[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:

enter image description here

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
$\endgroup$
10
  • $\begingroup$ "The response part of your edit should've been a comment". Nah, it was too long for a comment, I figured you'd see it anyway, or be notified of a change in an answer. $\endgroup$ – Cedron Dawg Sep 9 '20 at 15:11
  • $\begingroup$ Your continuous derivation appears invalid; tried it with a simple cosine - inverse yields negative of derivative in real part, but also an entire imaginary part (which should be zero), and normalization is also off. Problem seems to be in assuming $x[n] \rightarrow x(n)$; an inverse "inverts" the forward transformation, and a forward transform of a continuous signal would have $N=\infty$. $\endgroup$ – OverLordGoldDragon Sep 10 '20 at 13:03
  • $\begingroup$ How's it work with DTFT then? Not exactly sure, but I'd explain it as the encoding inserting as much information (i.e. $\omega$) 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. -- Lastly, it absurdly implies ability for infinite compression, as you're restoring an infinite $x(n)$ from $N$ finite coefficients. $\endgroup$ – OverLordGoldDragon Sep 10 '20 at 13:07
  • $\begingroup$ @OverLordGoldDragon This is a puzzler, I am looking at it. For the DTFT, you will have to center the spectrum on zero as -1 is not the same as N-1 in the interpolation. It helps that you used a degenerate case. One small point, I've added a ",False" parameter to your linspace call to not include the endpoint. Funny, since you mentioned that as an issue earlier. $\endgroup$ – Cedron Dawg Sep 10 '20 at 14:24
  • $\begingroup$ Good catch, and it is an issue; endpoints are akin to boundary conditions and shouldn't be treated lightly. The discrete derivative does accurately differentiate every case I've tried so far, but the slope addition seems to be a redundancy; DFT already assumes the input signal to be N-periodic, and where would we fetch $x[-1]$ anyway? The indices also should run from 1 to N - 1 instead of from 0, but I'm wondering how to 'justify' this restriction mathematically. $\endgroup$ – OverLordGoldDragon Sep 10 '20 at 14:34
2
$\begingroup$

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.

$\endgroup$
16
  • $\begingroup$ This is great - we're almost doing research here with this 'debate'. Maybe I should dedicate a Q&A to it, though, as we're going quite a bit off topic. $\endgroup$ – OverLordGoldDragon Sep 11 '20 at 10:41
  • $\begingroup$ @OverLordGoldDragon Actually I consider it a huge waste of my time. However, if it can drive a stake in the heart of the incredibly misleading assertion that "The DFT assumes the signal is periodic" for at least one person, it was worth it. $$ $$ I hope you can see that RB-J's "challenge" is pertinent to your question. The value of $X[k]\left( 1 - e^{-i\frac{2\pi}{N}k} \right)$ can be seen as the DFT of the original signal minus the original signal rotated by one sample. $\endgroup$ – Cedron Dawg Sep 11 '20 at 14:30
  • $\begingroup$ Can relate with your disposition, but mind that we're dealing with mathematics and there's no room for subjective disagreement; everything can be resolved exactly. 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. Explaining this, and explicitly clarifying when periodicity does and does not apply, may create a world-first complete primer on this topic (at least the freely-accessible one). $\endgroup$ – OverLordGoldDragon Sep 11 '20 at 14:37
  • $\begingroup$ Puzzle not solved; try x = np.random.randn(N). I'm unsure you followed my argument from information standpoint; your derivation implies absurdities. It's not a math typo, it's literally impossible. $\endgroup$ – OverLordGoldDragon Sep 11 '20 at 15:43
  • $\begingroup$ @OverLordGoldDragon That's mighty praise, but I doubt it deserves it. Two points to emphasize though: $$ $$ 1. Definitions are important $$ $$ 2. If you can't derive an assertion from definitions then you don't really understand it $$ $$ I'd also say there is room for subjectivity in Mathematics in the form of unproven conjectures and any results derived by assuming their truth (or falsehood). I would use the term "exhibit periodicity" rather than "demand periodicity". $\endgroup$ – Cedron Dawg Sep 11 '20 at 15:44
2
$\begingroup$

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.

$\endgroup$
13
  • $\begingroup$ I'm afraid you are a bit incorrect here as well. First, it would be nicer if you used $Y[k]$ for your differential as Olli and I did to distinquish it from the original DFT. Second, you are still missing the proper factor in your sine argument. $\sin(k/2)$ would have a ton of oscillations across your range. $\sin\left( \frac{\pi}{N}k \right)$ has the argument going from 0 to $\pi$ as $k$ goes from 0 to $N$. Likewise the arg will range from $ -\frac{\pi}{2}$ to $ \frac{\pi}{2}$ without any wrapping. $\endgroup$ – Cedron Dawg Sep 9 '20 at 14:41
  • $\begingroup$ @CedronDawg Right, I had the k/2 fixed - didn't actually derive it, but tossed |1 - e^(-ix)| to WolframAlpha and forgot about everything else. As for $Y[k]$, I am symbolically showing the direct effect on the non-differentiated $X[k]$. Also good lead on the alternative derivative; I'm going by Python's Numpy's implementation, which may not necessarily be best - but for this Q&A we rather finish the started. $\endgroup$ – OverLordGoldDragon Sep 9 '20 at 14:45
  • $\begingroup$ @CedronDawg Also, $k$ goes from $0$ to $N-1$ in one period, and the $\angle (0 + 0i)$ case needs accounting for. $\endgroup$ – OverLordGoldDragon Sep 9 '20 at 14:46
  • $\begingroup$ When talking about the general behavior, whether you include N or not isn't a big deal. So, yeah, you should have one lobe of the sine function, starting at 0, almost reaching 0 again at the end at N-1. Not two lobes as you show. Similar with the angle case, indeed arg(0) is undefined, but I think you will find a limit argument will quell any discomfort. Just put an open circle in the zero position and you'll be good. $\endgroup$ – Cedron Dawg Sep 9 '20 at 15:16
  • $\begingroup$ I didn't notice your scale. The graph is correct and I did my first comment incorrectly. The arg in my equation goes from $ -\frac{\pi}{2}$ to $ -3 \frac{\pi}{2}$. You can add $\pi$ (instead of $2\pi$) and still be correct because of the absolute value on the magnitude part. So, it goes from $ \frac{\pi}{2}$ to $ - \frac{\pi}{2}$ as shown in your graph. $\endgroup$ – Cedron Dawg Sep 9 '20 at 15:31
0
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ I'm a little bit disappointed in your implementation of my formula. Check out the latest followup in my first answer. ;-) $\endgroup$ – Cedron Dawg Sep 13 '20 at 12:18
  • $\begingroup$ Edited your code a little. So tell me, this makes sense to you? Because that's the consequence of your trick to get the nice numbers. $\endgroup$ – OverLordGoldDragon Sep 13 '20 at 13:38
  • 1
    $\begingroup$ @OlliNiemitalo Interesting, seem to meet all properties - I'll have a closer look later, thanks for the effort. $\endgroup$ – OverLordGoldDragon Sep 14 '20 at 16:45

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.