Forums

FFTW + low pass filtering + reverse FFTW

Started by greg838 April 6, 2006
Dear FFTW experts,

I'm using 2d r2c fftw on images, and now I would like to do a low pass
filtering. I compute the module
(sqrt(realpart*realpart+imagpart*imagpart)) of my data, and if it exceeds
a threshold, I set it to zero. Therefore, I get an array of real data, and
I don't know if I can use c2r transform by fulfilling my fftw_complex array
with my data (e. g. I put the value in the real part, an I set the
imaginary part to zero, but i feared i loose some datas), or if I must use
a r2r transform (which one, in this case)?

Sorry for my bad english.

Greetings
greg



greg838 skrev:
> Dear FFTW experts,
Have used the FFT, not the FFTW oas such...
> I'm using 2d r2c fftw on images, and now I would like to do a low pass > filtering. I compute the module > (sqrt(realpart*realpart+imagpart*imagpart)) of my data, and if it exceeds > a threshold, I set it to zero.
Is this what is called a "top hat filter"?
> Therefore, I get an array of real data, and > I don't know if I can use c2r transform by fulfilling my fftw_complex array > with my data (e. g. I put the value in the real part, an I set the > imaginary part to zero, but i feared i loose some datas), or if I must use > a r2r transform (which one, in this case)?
The main problem here is to maintain certain symmetries. First of all, you need to know where the DC component is. In some FFT implementations, it is in a corner, in others it is near the center. Assuming it is near the center, you need to have an as symmetric circle around DC as possible. If you are able to do that, the imaginary part of your IDFT'ed image ought to be negligable. Basically, this is a problem of book-keeping. If you can manage that, use the c2c IFFT and check the imaginary part. The magnitudes should be on the order of 1e-12 -- 1e-15, assuming you use double precision floating point numbers. Once you are confident you zero the correct numbers, you might save some run-time by using the c2r transform. But don't switch until you are confident you get the correct result! Rune
Ok I think I have succeed to do the stuff. I give my algorithm (in C++) if
other people are interested (and if I have made a mistake correct me, but
looking the result I think it seems to be quite right)

// m_out_p is my fft_complex array which contains the fftw_r2c results
//rad is the parameter of the low pass filter, the frequency response the
//frequency filter is : (k, l) are the horizontal and vertical frequency
//H(k,l) = 1 if sqrt(k&#2013266098;+l&#2013266098;) < threshold
//       = 0 if sqrt(k&#2013266098;+l&#2013266098;) >= threshold
//in fact,  rad equals threshold&#2013266098;

int ny = (int)floor( (double)ImageWidth/2) + 1; // Get the width of the
fftw_r2c array

	for(int j = 0; j < ny; j++)
	{
		for(int i = 0; i < ImageHeight; i++)
		{
			if ( rad < i*i + j*j )
			{
// Set to 0 real part and imaginary part
//in this test, i and j represent the horizontal and vertical frequency
				(*(m_out_p + (i*ny+j)))[0] = (double)0;
				(*(m_out_p + (i*ny+j)))[1] = (double)0;
			}

		}
	}

and fftw reverse.

Greetings,
greg
You can use the c2r transform, because *any* real-valued function of
the modulus of an r2c spectrum maintains the conjugate symmetry
required by complex-to-real (c2r) DFTs.

Regards,
Steven G. Johnson

in article 1144869202.536422.116380@t31g2000cwb.googlegroups.com,
stevenj@alum.mit.edu at stevenj@alum.mit.edu wrote on 04/12/2006 15:13:

> You can use the c2r transform, because *any* real-valued function of > the modulus of an r2c spectrum maintains the conjugate symmetry > required by complex-to-real (c2r) DFTs.
just curious, Steve, does the FFTW r2c (and c2r) put (and expect) the Nyquist compenent (the Re{X[N/2]}) into the Im{X[0]} bin? or does it do something else with it? -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
robert bristow-johnson wrote:
> just curious, Steve, does the FFTW r2c (and c2r) put (and expect) the > Nyquist compenent (the Re{X[N/2]}) into the Im{X[0]} bin? or does it do > something else with it?
No, it stores them as two complex numbers (with zero imaginary parts) at X[0] and X[N/2]. i.e. the r2c transform of a real array of length N is a complex array of length N/2+1. This wastes a miniscule amount of space, of course, but it generalizes much better to higher-dimensional transforms. (If you try to use packing tricks with higher-dimensional r2c FFTs you get absurdly complicated data formats...it's not worth it in my opinion.) We have an alternate 1d FFT interface that uses a different format which doesn't have any wasted space: the real parts followed by the imaginary parts in reverse order (the FFTW_R2HC transform). This is not really to save space, however; it is because the real-input-specialized Cooley-Tukey algorithm naturally works nicely with this format (which was first proposed by Sorensen in the 80's if I recall correctly). FFTW uses this format internally for real-input transforms of odd N size. Cordially, Steven G. Johnson
Thank you all for your replies.

Greetings,
greg