Calculating the first derivative of an image using DFT
up vote
2
down vote
favorite
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.
dft python
asked Nov 13 at 22:31
RonaldB
1122
add a comment |
up vote
2
down vote
favorite
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.
dft python
asked Nov 13 at 22:31
RonaldB
1122
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.
– Andy Walls
Nov 13 at 22:51
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/…
– Andy Walls
Nov 13 at 23:08
add a comment |
up vote
2
down vote
favorite
up vote
2
down vote
favorite
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.
dft python
asked Nov 13 at 22:31
RonaldB
1122
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.
dft python
dft python
asked Nov 13 at 22:31
RonaldB
1122
asked Nov 13 at 22:31
RonaldB
1122
asked Nov 13 at 22:31
RonaldB
1122
asked Nov 13 at 22:31
RonaldB
1122
asked Nov 13 at 22:31
RonaldB
1122
1122
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.
– Andy Walls
Nov 13 at 22:51
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/…
– Andy Walls
Nov 13 at 23:08
add a comment |
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.
– Andy Walls
Nov 13 at 22:51
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/…
– Andy Walls
Nov 13 at 23:08
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.
– Andy Walls
Nov 13 at 22:51
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.
– Andy Walls
Nov 13 at 22:51
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/…
– Andy Walls
Nov 13 at 23:08
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/…
– Andy Walls
Nov 13 at 23:08
add a comment |
1 Answer
1
active
oldest
votes
up vote
3
down vote
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^{jomega}) - e^{-j omega} X(e^{jomega}) =X(e^{jomega})(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 jOmega 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.
answered Nov 13 at 23:08
Fat32
14k31128
add a comment |
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
up vote
3
down vote
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^{jomega}) - e^{-j omega} X(e^{jomega}) =X(e^{jomega})(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 jOmega 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.
answered Nov 13 at 23:08
Fat32
14k31128
add a comment |
up vote
3
down vote
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^{jomega}) - e^{-j omega} X(e^{jomega}) =X(e^{jomega})(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 jOmega 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.
answered Nov 13 at 23:08
Fat32
14k31128
add a comment |
up vote
3
down vote
up vote
3
down vote
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^{jomega}) - e^{-j omega} X(e^{jomega}) =X(e^{jomega})(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 jOmega 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.
answered Nov 13 at 23:08
Fat32
14k31128
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^{jomega}) - e^{-j omega} X(e^{jomega}) =X(e^{jomega})(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 jOmega 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.
answered Nov 13 at 23:08
Fat32
14k31128
edited Nov 13 at 23:18
answered Nov 13 at 23:08
Fat32
14k31128
answered Nov 13 at 23:08
Fat32
14k31128
answered Nov 13 at 23:08
Fat32
14k31128
14k31128
add a comment |
add a comment |
Thanks for contributing an answer to Signal Processing Stack Exchange!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
Some of your past answers have not been well-received, and you're in danger of being blocked from answering.
Please pay close attention to the following guidance:
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
To learn more, see our tips on writing great answers.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fdsp.stackexchange.com%2fquestions%2f53290%2fcalculating-the-first-derivative-of-an-image-using-dft%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
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.
– Andy Walls
Nov 13 at 22:51
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/…
– Andy Walls
Nov 13 at 23:08