2
$\begingroup$

I need to calculate the first derivative of a greyscale image (a 2D array) using a DFT function I built (which works). Unfortunately, the results don't seem to be correct.

In the fourier domain, the derivative d/dx is given as F(u,v)* 2*pi*i/N * u, where u is the x-axis transformed, N is the size of one of the matrix's dimensions, v being the other one.

Attached is the code. What bothers me is that I'm not getting the same results as I would differentiating by convolution with (1,-1) or (1,-1) as a column vector.

def derivative(fourier_signal):
"""
Derivative in fourier domain is multiplying by u or v, and 2pi*i/N
:param fourier_signal:
:return:
"""
N = np.shape(fourier_signal)[ZERO]
M = np.shape(fourier_signal)[ONE]
u = np.arange(N)
v = np.arange(M)
du = fourier_signal * (u*TWO_PI*1j)/N
dv = fourier_signal * (v*TWO_PI*1j)/M
return du, dv


def fourier_der(im):
    # Calculate DFT2
    dft_image = DFT2(im)
    # Function that Multiply by rows by u, columns by y
    du, dv = derivative(dft_image)
    shifted_du, shifted_dv = np.fft.fftshift(du), np.fft.fftshift(dv)
    dx, dy = IDFT2(shifted_du), IDFT2(shifted_dv)

I'm not looking for easy answers on how to do it, but rather a direction to as to why my output is incorrect.

$\endgroup$
2
  • $\begingroup$ 1) Aliasing. The ideal derivative filter you're using in the frequency domain is infinitely long in the time domain. You're getting aliases in the time domain. 2) [1, -1] as filter taps is a scaled truncation of the ideal derivative filter. It is in error at high frequencies. So in light of 1 & 2, you should really specify over which frequencies how much error you can tolerate. Then you can build a filter to that specification. $\endgroup$
    – Andy Walls
    Nov 13, 2018 at 22:51
  • $\begingroup$ BTW, GNURadio has a little stand alone utility for building Minimum Mean Squared Error differentiating FIR filter taps. GNURadio uses it to build a polyphase filter bank for an interpolating differentiator, with 8 taps in each polyphase filter arm, and MMSE over the frequencies in the interval $[-f_s/4, f_s/4]$. github.com/gnuradio/gnuradio/tree/master/gr-filter/lib/… $\endgroup$
    – Andy Walls
    Nov 13, 2018 at 23:08

1 Answer 1

3
$\begingroup$

They are not the same. Using 1D notation,

the discrete-time (backward) first difference is $x[n] - x[n-1]$ whose frequency domain DTFT equivalent is

$$ x[n]-x[n-1] \leftrightarrow X(e^{j\omega}) - e^{-j \omega} X(e^{j\omega}) =X(e^{j\omega})(1- e^{-j \omega}) $$

which becomes $$ x[n]-x[n-1] \longleftrightarrow X[k](1 - e^{-j \frac{2 \pi}{N} k})$$ using the DFT to implement it.

The FIR impulse response of the discrete-time system that implements the first difference is therefore,

$$ h[n] = \delta[n] - \delta[n-1]$$

The first derivative of a continuous-variable function $x(t)$ is $x'( t)$ and in CTFT domain it becomes : $$ x'(t) \longleftrightarrow j\Omega X(\Omega) $$

where the analog system frequency response is

$$H_d(\Omega) = j \Omega $$

which is not implementable in digital form, but a bandlimited approximation to it is attained under a sampling period of $T$ that yields an equivalent discrete-time frequency response of a (bandlimited) differentiator as

$$ H_d(e^{j \omega}) = j \frac{\omega}{T} ~~~, ~~~~ |\omega| < \pi$$

The associated discrete-time (IIR) impulse response is $$ h_d[n] = \text{IDTFT} \{ H_d(e^{j \omega}) \} = \text{IDTFT} \{ j \frac{\omega}{T} \} $$

Practically you would truncate and window $h_d[n]$ before using.

So they are not the same.

$\endgroup$

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.