Low-pass filtering via FFTW: artefacts and wrong scaling

Dear members,
I work on magneto-optical imaging (MOIF) of magnetic nanostructures. To
compare exact calculations of magnetic field distributions with
experimental data, I want to perform some low-pass filtering on my
theoretical simulations. As a first try I convolved my data with a
gaussian kernel via complex multiplication in fourier domain using 2D
complex-to-complex dft (FFTW release 3.1.1 included in C++ code; g++
compiler on Unix system).

My data is represented as a greyscale image (grey value corresponds to
magnetic field strength), e.g. a bright (i.e magnetized) quadratic pattern
in a dark (i.e. nonmagnetic) surrounding. For the convolution I centered
such a quadratic pattern and a normalized gaussian bell each on a
1000x1000 point grid, performed the forward trafos, complex multiplication
and subsequent back trafo.

The output shows the favoured features like an unsharpening of the
pattern's edges, but two problems occur:

- the original pattern is cut in four congruent pieces, each of them
displaced to one corner of the underlying grid

- the outcoming field strengths, i.e. the intensity of the greyscale image
is several orders of magnitude too large, even though I scaled the output
with 1/(N*N) (where N=1000 is the number of grid points per axis). I feel
like this point has something to do with the chosen variance of the
normalized gaussian.

I'd be very thankful for any help with this.
Greetings, J.Norpoth

"norpoth" <norpoth@ump.gwdg.de> wrote in message 
news:v7SdnenlaPl1T83ZnZ2dnUVZ_sKdnZ2d@giganews.com...
> Dear members,
> I work on magneto-optical imaging (MOIF) of magnetic nanostructures. To
> compare exact calculations of magnetic field distributions with
> experimental data, I want to perform some low-pass filtering on my
> theoretical simulations. As a first try I convolved my data with a
> gaussian kernel via complex multiplication in fourier domain using 2D
> complex-to-complex dft (FFTW release 3.1.1 included in C++ code; g++
> compiler on Unix system).
>
> My data is represented as a greyscale image (grey value corresponds to
> magnetic field strength), e.g. a bright (i.e magnetized) quadratic pattern
> in a dark (i.e. nonmagnetic) surrounding. For the convolution I centered
> such a quadratic pattern and a normalized gaussian bell each on a
> 1000x1000 point grid, performed the forward trafos, complex multiplication
> and subsequent back trafo.
>
> The output shows the favoured features like an unsharpening of the
> pattern's edges, but two problems occur:
>
> - the original pattern is cut in four congruent pieces, each of them
> displaced to one corner of the underlying grid
>
> - the outcoming field strengths, i.e. the intensity of the greyscale image
> is several orders of magnitude too large, even though I scaled the output
> with 1/(N*N) (where N=1000 is the number of grid points per axis). I feel
> like this point has something to do with the chosen variance of the
> normalized gaussian.

Maybe this thread will give you some hints:

http://groups.google.com/group/comp.dsp/browse_thread/thread/e3c0bcc5bd09ae7c/3a37e015d004f554?lnk=st&q=comp.dsp+%22Fred+Marshall%22+gaussian+image&rnum=3&hl=en#3a37e015d004f554 

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message 
news:-e6dnfGoRbkhQM3ZRVn-sg@centurytel.net...

>
> Maybe this thread will give you some hints:
>
> http://groups.google.com/group/comp.dsp/browse_thread/thread/e3c0bcc5bd09ae7c/3a37e015d004f554?lnk=st&q=comp.dsp+%22Fred+Marshall%22+gaussian+image&rnum=3&hl=en#3a37e015d004f554
>

I might have also said:

I'm going to assume that nonvarying energy is in the middle of the grid 
(around 500,500).
For me, it's just that much easier to deal with the information - to read it 
and visualize it, to plot it, etc.

This requires that you do a "shift" in the frequency domain:

See:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fftshift.html

Then, you apply a lowpass filter with a peak, the passband, at the *center* 
of the frequency array.

Where scaling is concerned you may want the Gaussian in frequency to be 
normalized so that it has a peak of 1.0 at the center - thus not changing 
the energy of components at low frequencies.

Then you shift the information back before the ifft.

Or, it may be more efficient to do this:

Get the filter you want, centered at 500,500.  Then shift the filter and use 
its shifted version on non-shifted data ffts.  You can work out the details 
and test it.

Then, with FFTW, you have to do all the scaling - as it does none I believe.

This could mean scaling by 1/N^2 with the ifft. Or it could mean scaling by 
1/N in the forward fft and in the ifft.  NOTE: these apply to 2D and not 1D!

Fred 

Thanks for the quick answer. Since we have a public holiday here in Germany
I won't give it a go until next week. I'll report on progress.
Greetings, J.Norpoth

Hello,
the first problem is solved. After performing the pointwise complex
multiplication of the two convolution partner in fourier domain I had to
multiply each of these products by (-1)^(i+j), where i,j label the
respective position on the underlying NxN grid. According to the FFTW FAQ
(Question 3.10) this trick has to be applied to the input of a forward
trafo to shift the zero-frequency component to the center of the output
and it seems to work in the other direction too.
Unfortunately the scaling problem remains.

Greetings, J.Norpoth

Hi,

Was the scaling problem resolved?

Regards