Disclaimer - Edit September 16, 2015
This message has changed significantly due to a small but fundamental flaw in my understanding of the derivative operation in the frequency domain. Most of this post has changed since the previous iteration.
I would like to thank Luis Mendo for helping me debug why this did not work on other images besides the one that was presented in the original version of the message.
When it comes to discrete time data, there is no such thing as a derivative. A derivative is an analytical tool that exists only in continuous time . In discrete time, we can only navigate through differences . So I'm going to suggest that you mean a difference operation, not a derivative. One of the most common approximations to a derivative is to use the direct difference operation. Specifically, for a 1D discrete sequence x[n] output sequence y[n] , where the direct difference operation is applied, is defined as:
y[n] = x[n+1] - x[n], for n = 1, 2, ..., M-1
M is the total number of samples in your discrete sequence. Naturally, the value will be less due to the nature of the forward diversity operation.
To do the same in the frequency domain, you must use the shift version from the properties of the discrete Fourier transform. Specifically, given the version of the DFTed signal x[n] , which is X(k) , the following property connects the time domain and the frequency domain:
Here i is the complex number or square root of -1, and k is the angular frequency in the frequency domain. F^{-1} is the inverse Fourier transform of the signal X(k) , which is the Fourier transform of the signal in the time domain x[n] . M is also the length of the signal x[n] , and M is the amount of shift by which you want to move the signal. Therefore, this means that if you want to calculate the shift operation by shifting the positions of M , you just need to accept the Fourier transform, multiply each component by exp(-i*2*pi*m*k/M) , where k is the index to the Fourier transform point and take the inverse Fourier transform of this intermediate result.
We emphasize that the last element of the Fourier transform is meaningless, because the difference operation leads to a signal with one smaller element, so it is important to remove the last element from this result.
In view of the foregoing, in order to calculate the approximation to the derivative, we need to find the Fourier transform y[n] or Y(k) and can be calculated as follows:
Note that the shift is -1 such that x[n+1] = x[n-(-1)] . Therefore, you must calculate the Fourier transform, multiply each corresponding value by exp(i*2*pi*k/M) - 1 and take the inverse.
It is not much to go into 2D. Here I assume that we are performing a difference operation in both directions - horizontal and vertical. Remember that we have two variables that represent horizontal and vertical spatial coordinates. We also call this the spatial domain, not the time domain, since the variables correspond to our position in the 2D space. Following the wording above, switching to 2D is very simple:
Here u and v represent the spatial coordinates of the 2D-discrete difference operation y[u,v] . u and v are two-dimensional frequency components, and M and N are the number of columns and rows of a two-dimensional signal. Pay particular attention to the fact that you are actually performing a horizontal difference operation followed by a vertical difference operation. In particular, you perform the horizontal difference operation for each row independently, followed by the vertical difference operation for each column separately. Actually, the order doesn't matter. You can do one after the other, in any order you choose. It is important to note that the 2D Fourier transform of the signal remains unchanged, and you simply multiply each value by some complex constant. However, you need to make sure that you delete the last row and last column of the result, because each operation results in a signal that is one point less than the original length of each row and each column due to the observed difference observation that we talked about,
Assuming you already have FFT functions, you just need to take each two-dimensional spatial coordinate and multiply by (exp(i*2*pi*U/M) - 1)*(exp(i*2*pi*V/N) - 1) . For each 2D FFT component stored in (U,V) , we use the same frequency coordinates and multiply this location by the product of these two things before. After that, you should take the inverse FFT, and we compare this with the actual derivative of the spatial domain. We must see that they are the same.
Note that when performing 2D FFT, the frequencies are normalized , so that they span between [-1,1] in both dimensions. We must take this into account when mapping the equivalence of derivative operations in both domains. In addition, to make things a lot easier, it would also be wise if you would shift the frequency components so that they appear in the center of the image. When performing 2D FFT, the components are positioned so that the DC component or source is located in the upper left corner of the image. Move everything to the center using fftshift .
Now let's start with something small. For this procedure, we assume that the size of each measurement is odd in order to simplify and unambiguously change the frequency offset through fftshift . Say we had a 101 × 101 gauss filter with a standard deviation of 10, and we take the FFT of that, then fftshift . So:
h = fspecial('gaussian', [101 101], 10); hf = fftshift(fft2(h));
If we want to find the derivative, we define a grid of frequency points that span the image in the frequency domain, starting with 0 as the source, which is now the center . We can do this through a meshgrid , and so:
[U,V] = meshgrid(-50:50, -50:50);
In general, these horizontal and vertical coordinates span from [-floor(cols/2),floor(cols/2)] for horizontal and [-floor(rows/2),floor(rows/2)] for vertical. Please note that the above grid has the same dimensions as the image, but the center is at (0,0) and we occupy from -50 to 50 in both dimensions.
Now perform a difference operation in the frequency domain in both directions. Let's do it for both x and y , so the first derivative in both directions:
hf2 = (exp(1i*2*pi*U/101) - 1).*(exp(1i*2*pi*V/101) - 1).*hf;
Note that we had to normalize the frequencies by dividing by 101 due to the fact that the FFT frequencies are normalized. 101 x 101 due to image size. If we want this back in the time domain, just take the inverse FFT. However, before accepting the inverse FFT, we need to undo the shift we made with ifftshift . We also use real to eliminate any of the residual complex-valued components due to numerical errors. You will see that the result is comprehensively evaluated, but the imaginary components are so small in magnitude that we can ignore them. Thus:
out_freq = real(ifft2(ifftshift(hf2)));
Remember to also remove the last row and column:
out_freq = out_freq(1:end-1,1:end-1);
If you want to compare this output with the time difference operation of the filter h , we can do this with diff over the rows, then using this result, we go through the columns. Thus:
out_spatial = diff(diff(h, 1, 1), 1, 2);
The second parameter indicates the order in which the difference operation operates. This is a first-order difference, so the horizontal and vertical difference operations are set to 1. The third parameter indicates in which dimension you are making the difference. First we do this column by column, and then apply the intermediate result row by row. Then we can show how they look:
figure; subplot(1,2,1); imshow(out_spatial, []); title('Derivative from spatial domain'); subplot(1,2,2); imshow(out_freq, []); title('Derivative from frequency domain');
And we get:
Cool! This is consistent with what I know about the Gaussian derivative.
Bonus
As a bit of a bonus, I would like to point you to some great slides from Cornell University: http://www.cs.cornell.edu/courses/cs6670/2011sp/lectures/lec02_filter.pdf
In particular, this is what looks Gaussian if you want to look at it in 3D, where we have both horizontal and vertical coordinates, and the Z value is the height of the Gaussian in 3D:
So, this is what happens if we take the derivative of both directions independently. The top shows what happens spatially, and the bottom shows you what will happen if we look at it directly above the bird's eye:
The most negative part of the surface is drawn in black, and the most positive part of the surface is drawn in white, and everything else is scaled in intensity between them ... so if we take the spatial derivatives in one direction, then use this intermediate result and take the spatial derivative in the other direction , you will see that we get what we saw above. Play with different M and N values and see what you get. This property, which I have not seen, is often used in practice, but it is certainly useful to know.