Calculating the first derivative of an image using DFT
$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.
dft python
asked Nov 13 '18 at 22:31
RonaldBRonaldB
1122
$endgroup$
add a comment |
$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.
dft python
asked Nov 13 '18 at 22:31
RonaldBRonaldB
1122
$endgroup$
$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 '18 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 '18 at 23:08
add a comment |
$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.
dft python
asked Nov 13 '18 at 22:31
RonaldBRonaldB
1122
$endgroup$
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 '18 at 22:31
RonaldBRonaldB
1122
asked Nov 13 '18 at 22:31
RonaldBRonaldB
1122
asked Nov 13 '18 at 22:31
RonaldBRonaldB
1122
asked Nov 13 '18 at 22:31
RonaldBRonaldB
1122
asked Nov 13 '18 at 22:31
RonaldBRonaldB
1122
1122
$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 '18 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 '18 at 23:08
add a comment |
$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 '18 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 '18 at 23:08
$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 '18 at 22:51
$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 '18 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 '18 at 23:08
$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 '18 at 23:08
add a comment |
1 Answer
1
active
oldest
votes
$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^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 frac2 piN 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 fracomegaT ~~~, ~~~~ |omega| < pi$$
The associated discrete-time (IIR) impulse response is
$$ h_d[n] = textIDTFT H_d(e^j omega) = textIDTFT j fracomegaT $$
Practically you would truncate and window $h_d[n]$ before using.
So they are not the same.
answered Nov 13 '18 at 23:08
Fat32Fat32
15.1k31231
$endgroup$
add a comment |
Your Answer
StackExchange.ifUsing("editor", function ()
return StackExchange.using("mathjaxEditing", function ()
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix)
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
);
);
, "mathjax-editing");
StackExchange.ready(function()
var channelOptions =
tags: "".split(" "),
id: "295"
;
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function()
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled)
StackExchange.using("snippets", function()
createEditor();
);
else
createEditor();
);
function createEditor()
StackExchange.prepareEditor(
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader:
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
,
noCode: true, onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
);
);
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
1 Answer
1
active
oldest
votes
1 Answer
1
active
oldest
votes
active
oldest
votes
active
oldest
votes
$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^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 frac2 piN 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 fracomegaT ~~~, ~~~~ |omega| < pi$$
The associated discrete-time (IIR) impulse response is
$$ h_d[n] = textIDTFT H_d(e^j omega) = textIDTFT j fracomegaT $$
Practically you would truncate and window $h_d[n]$ before using.
So they are not the same.
answered Nov 13 '18 at 23:08
Fat32Fat32
15.1k31231
$endgroup$
add a comment |
$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^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 frac2 piN 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 fracomegaT ~~~, ~~~~ |omega| < pi$$
The associated discrete-time (IIR) impulse response is
$$ h_d[n] = textIDTFT H_d(e^j omega) = textIDTFT j fracomegaT $$
Practically you would truncate and window $h_d[n]$ before using.
So they are not the same.
answered Nov 13 '18 at 23:08
Fat32Fat32
15.1k31231
$endgroup$
add a comment |
$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^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 frac2 piN 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 fracomegaT ~~~, ~~~~ |omega| < pi$$
The associated discrete-time (IIR) impulse response is
$$ h_d[n] = textIDTFT H_d(e^j omega) = textIDTFT j fracomegaT $$
Practically you would truncate and window $h_d[n]$ before using.
So they are not the same.
answered Nov 13 '18 at 23:08
Fat32Fat32
15.1k31231
$endgroup$
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 frac2 piN 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 fracomegaT ~~~, ~~~~ |omega| < pi$$
The associated discrete-time (IIR) impulse response is
$$ h_d[n] = textIDTFT H_d(e^j omega) = textIDTFT j fracomegaT $$
Practically you would truncate and window $h_d[n]$ before using.
So they are not the same.
answered Nov 13 '18 at 23:08
Fat32Fat32
15.1k31231
edited Nov 13 '18 at 23:18
answered Nov 13 '18 at 23:08
Fat32Fat32
15.1k31231
answered Nov 13 '18 at 23:08
Fat32Fat32
15.1k31231
answered Nov 13 '18 at 23:08
Fat32Fat32
15.1k31231
15.1k31231
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.
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
$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 '18 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 '18 at 23:08