Finding equivalent gaussian filter mask in frequency domain for given mask in space domain

3
votes

So far I've implemented a gaussian blur filter entirely in the space domain, making use of the separability of the gaussian, that is, applying a 1D gaussian kernel along the rows and then along the columns of an image. That worked fine.

Now, given only with the size N of the NxN convolution matrix of the space domain, I want to achieve the exact same blurred image over the frequency domain. That means that I'll load the image into a matrix (numpy, I'm using python), apply the FFT on it (then I have G(x,y)), and then I have to have also a filter H(u,v) in the frequency domain that also resembles the shape of some 2d gaussian, with its center value being 1.0 and then having values falling off to 0 the further away from the center I am. I do then the multiplication in frequency domain (before I have to consider to do a center-shift of H) and then apply the iFFT.

The trouble I have is to find the exact formula (i.e. to find sigma, the std-deviation) that will result in the corresponding H(u,v). From the space domain, if I have been given a mask-size N, I know that the std-dev sigma can be approximated as sigma=(maskSize-1)/2/2.575, e.g. for a mask size N=15 I get std-dev=2.71845 for e^-(x²/2sigma²), just considering 1D cases for now.

But how do I get sigma for the frequency domain?

Funny thing is btw that in theory I know how to get sigma, using Mathematica, but the result is pure bogus, as I can demonstrate here:

gauss1d[x_, sigma_] := Exp[-(x^2)/(2 sigma^2)]
Simplify[FourierTransform[gauss1d[x, sigma], x, omega], sigma > 0]

The result is E^(-(1/2) omega^2 sigma^2) * sigma

This is bogus because it turns, in the exponent of the E function, the 1/sigma² into a sigma². Consequently, if you draw this, you will see that the standard deviation has become a lot smaller, since the H(u,v)-gaussian is a lot "thinner". However, it should actually be a lot wider in the frequency domain than in the space domain!! It doesn't make any sense...

1
mathfftblurfrequencygaussian
Why can't you just take the FFT of the spatial domain Gaussian filter and apply that to your frequency domain image ? – Paul R
I tried that. I created the spatial domain gaussian filter, whose matrix has the same dimensions as the image I want to blur. Then, after applying the FFT of it, I'll unfortunately get "ripple effects" in the result. That means, I do get the FFT, but it also has riples along the x and y axis (not, however, along the diagonals). Of course I could manually process the array to remove those, but that looks like a hack to me. I think the ripples are an artifact due to discrete FFT computation. – NameZero912
See here: i56.tinypic.com/21d7o0.jpg – NameZero912
Ripples are reasonable. Gaussian kernel has infinite domain, and you restricted it to certain box. That's effectively convolving with the boxcar's fourier image, which has many ripples. – mbaitoff
There's nothing wrong with your Fourier transform of the Gaussian function. Try reading the Wikipedia article on the FT and it should prove enlightening. Also, you forgot to scale the spatial Gaussian by 1 / sqrt(2*pi*sigma^2) – Mike Bailey

1 Answers

8
votes

The Fourier Transform of the a Gaussian is a Gaussian as you can see from

http://en.wikipedia.org/wiki/Fourier_transform

But note that the std. dev. DOES invert!!!!

You say that is bogus BUT it is not. The frequency domain is, in some sense, the inversion of the time domain.

freq = 1/time

The standard deviation given is in time, when you transform it is still in time(the constant does not get transformed).

Suppose your you found the time version of the Gaussian using some s in terms of time. You transform the data into freq space. You can use that s and it will behave exactly the way it is suppose to. e.g., if you have a small s then it will cause the freq std. dev. in the frequency version to be large.

Again, this is because frequency is the inversion of time(again, in a sense).

Suppose your Gaussian has very small std. dev. Then it approximates a dirac delta function. We know this because gets transformed into a sinusoid in the freq domain. i.e., something that spans the whole frequency domain. (i.e., has infinite std. dev.(if it were a Gaussian).

Think of it like this: You are wanting to smooth in the frequency domain. Smooth what? High frequency components, right? By convoluting with a gaussian you are smoothing near by data. If the std. dev. is small you are keeping higher frequencies. In the frequency domain this means you are KEEPING more frequencies. But the conv is a multiplication in frequency domain. If we mulitplied by a thin gaussian in the frequency domain we would be left with a small group of frequencies.

G(t)*f(t) G[w]*f[w]

The first, a convolution. For a smooth filter we want G(t) to be "large"(std. dev. to be large). This means in we want less of the high frequency components(a sort of low pass filter). In the freq. domain we are multiplying by G[w]. This means that G[w] must be thin(and centered around the origin) so that we block out the highs.

I think basically you are not realizing that in the time domain we have a convolution and the frequency domain it is a multiplication. G cannot be the same in both. If G is thin in the time domain and thin in the frequency domain they will not result in the same effect. G thin in the convolution gives almost no effect but G thin in the freq. domain almost completely removes all the frequencies(just the very low are kept).