Forums

FFTW image processing frequency value question...(ad nauseum)

Started by dpjones0 September 16, 2008
Hi,

I have been using FFTW to do some image processing, including cross
correlations, etc.  I have found it to be very useful.  I've used the r2c/
c2r combination as well as forward and reverse fftw_plan_dft_1d, based on
examples I have found online.  I have been reading this discussion group
that has helped me figure a lot of things out.  However, in my latest
application, I still seem to be stuck on one issue.  

I really need to know the exact frequency values that each index in the
FFT resultant complex represents.  At this point, speed is not that much of
an issue, as much as understanding what I'm doing.  So for example, I set
the input data to an NxN (e.g. 64x64) size image and use an out-of-place
fftw_plan_dft_1d setting the input imaginary data to 0 and the real data to
the pixel values.  I then end up with an NxN size FFTW complex, where the
DC value is at index 0,0.  Does this also imply that in the unshifted FFTW
frequency domain data, the first entire column has a horizontal frequency
of 0 and the first entire row has vertical frequency of 0???.  I presumed
this was the case, since the index 0,0 corresponded to 0 vertical freq and
0 horizontal freq, but at this point it's not clear to me.  Does this mean
that for instance, along row=0, the horizontal frequency would go from:

(row,col)       horiz. freq.       vert. freq.
(0,0)             0                      0    \
(0,1)             1*1/N                  0     |-> 32 indices f=0 to
.5-1/N
(0,N/2 - 1)       (N/2-1)*1/N            0    /

(0,N/2)           - (N/2)*1/N            0\ 
(0,N-1)           - (1/(N/2)             0/  --> 32 indices f= -.5 to
-1/N

The reason I ask this is I got confused by trying to center the DC
component.  To look at power spectrums, etc with DC centered, in this
discussion group and elsewhere online, it is suggested to move the DC
component to N/2,N/2, and swap the quadrants (e.g. 1 and 4, 2 and 3).  

So, after doing that the DC component is now at row=N/2, col=N/2 (e.g.
32,32 for my 64x64 input image).  Then, the horizontal frequencies increase
from coordinate row=32,col=32 to coordinate row=32,col=63 essentially from
DC (0 freq) to PI/32 in 32 pixels.  However, the conjugate symmetric data
for the negative frequencies now occurs from coordinates row=31,col=31 to
row=31,col=1, with an extra value at row=31, col=0.  Why is the conjugate
symmetric data offset by a row?  What are the frequencies for each of the
points surrounding the center DC position row=32,col=32?

Alternatively, it was also posted that you could just multiply the input
image pixels at coordinates (i,j) by (-1)^(i+j) for instance.  The problem
that I have is that these two options do not present the same results.  If
I do that, there is no DC value at N/2,N/2.  In other words, the data from
row=32,col=32 to row=32,col=63 equals the conjugate symmetric data from
row=32,col=31 to row=32,col=1, with an oddball number in the 0 column.

I think I'm missing something fundamental, b/c this seems a lot more
confusing than it should be...


Thanks for any help in these matters.
-Dan
dpjones0


On Sep 16, 4:06&#2013266080;pm, "dpjones0" <dpjon...@comcast.net> wrote:
> Hi, > > I have been using FFTW to do some image processing, including cross > correlations, etc. &#2013266080;I have found it to be very useful. &#2013266080;I've used the r2c/ > c2r combination as well as forward and reverse fftw_plan_dft_1d, based on > examples I have found online. &#2013266080;I have been reading this discussion group > that has helped me figure a lot of things out. &#2013266080;However, in my latest > application, I still seem to be stuck on one issue. &#2013266080; > > I really need to know the exact frequency values that each index in the > FFT resultant complex represents. &#2013266080;At this point, speed is not that much of > an issue, as much as understanding what I'm doing. &#2013266080;So for example, I set > the input data to an NxN (e.g. 64x64) size image and use an out-of-place > fftw_plan_dft_1d setting the input imaginary data to 0 and the real data to > the pixel values. &#2013266080;I then end up with an NxN size FFTW complex, where the > DC value is at index 0,0. &#2013266080;Does this also imply that in the unshifted FFTW > frequency domain data, the first entire column has a horizontal frequency > of 0 and the first entire row has vertical frequency of 0???. &#2013266080;I presumed > this was the case, since the index 0,0 corresponded to 0 vertical freq and > 0 horizontal freq, but at this point it's not clear to me. &#2013266080;Does this mean > that for instance, along row=0, the horizontal frequency would go from: > > (row,col) &#2013266080; &#2013266080; &#2013266080; horiz. freq. &#2013266080; &#2013266080; &#2013266080; vert. freq. > (0,0) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 0 &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0 &#2013266080; &#2013266080;\ > (0,1) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 1*1/N &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0 &#2013266080; &#2013266080; |-> 32 indices f=0 to > .5-1/N > (0,N/2 - 1) &#2013266080; &#2013266080; &#2013266080; (N/2-1)*1/N &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0 &#2013266080; &#2013266080;/ > > (0,N/2) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; - (N/2)*1/N &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0\ > (0,N-1) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; - (1/(N/2) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 0/ &#2013266080;--> 32 indices f= -.5 to > -1/N > > The reason I ask this is I got confused by trying to center the DC > component. &#2013266080;To look at power spectrums, etc with DC centered, in this > discussion group and elsewhere online, it is suggested to move the DC > component to N/2,N/2, and swap the quadrants (e.g. 1 and 4, 2 and 3). &#2013266080; > > So, after doing that the DC component is now at row=N/2, col=N/2 (e.g. > 32,32 for my 64x64 input image). &#2013266080;Then, the horizontal frequencies increase > from coordinate row=32,col=32 to coordinate row=32,col=63 essentially from > DC (0 freq) to PI/32 in 32 pixels. &#2013266080;However, the conjugate symmetric data > for the negative frequencies now occurs from coordinates row=31,col=31 to > row=31,col=1, with an extra value at row=31, col=0. &#2013266080;Why is the conjugate > symmetric data offset by a row? &#2013266080;What are the frequencies for each of the > points surrounding the center DC position row=32,col=32? > > Alternatively, it was also posted that you could just multiply the input > image pixels at coordinates (i,j) by (-1)^(i+j) for instance. &#2013266080;The problem > that I have is that these two options do not present the same results. &#2013266080;If > I do that, there is no DC value at N/2,N/2. &#2013266080;In other words, the data from > row=32,col=32 to row=32,col=63 equals the conjugate symmetric data from > row=32,col=31 to row=32,col=1, with an oddball number in the 0 column. > > I think I'm missing something fundamental, b/c this seems a lot more > confusing than it should be... > > Thanks for any help in these matters. > -Dan > dpjones0
Consider an 8x8 example. Consider just the FFT of one row. You'll have: 0f, 1f, 2f, 3f, 4f, -3f, -2f, -1f You go from low frequency to high frequency and then back to low (considering -3f to be higher than -1f). All of the rows will have that frequency arrangement. Now consider the FFT of one column, You'll have: 0f 1f 2f 3f 4f -3f -2f -1f All the columns will have the same frequency arrangement. Now consider the 4,4 point in the output - that's the (4f,4f) location - your highest frequency. The 0,0 point is DC - your lowest frequency. The 3,3 point is (3f,3f), 2,2 is (2f,2f), etc. But because of the negative frequencies, the point 7,7 is actually (-1f,-1f). Similarly, point 6,7 is (-2f,-1f). Before centering the DC point, the 4,4 point is the highest frequency, and you go from high frequencies to low frequencies as you move towards the edges. When you center the DC output, you are in fact exchanging the data in the 4 quadrants. For N even, each one of those quadrants contains (N/ 2)x(N/2) points. For the N = 8x8 case, with the DC output centered (or rotated) to 4,4, the surrounding points are (row, column): (-1f, -3f) (0f, -1f) (1f, -1f) (-1f, 0f) DC (1f, 0f) (-1f, 1f) (0f, 1f) (1f, 1f) OW! My head hurts! It really is confusing. For instance, before centering, columns 0, 1 and 2 have the following frequency arrangements: 0f0f 1f0f 2f0f 0f1f 1f1f 2f1f 0f2f 1f2f 2f2f 0f3f 1f3f 2f3f 0f4f 1f4f 2f4f 0f-3f 1f-3f 2f-3f 0f-2f 1f-2f 2f-2f 0f-1f 1f-1f 2f-1f After centering, the arrangement gets a lot more confusing because of the quadrant swapping. I've also seen the multiply outputs bit, and I think it actually does work, because I saw some code that used it. But I've always just used 2 'for' loops to get it done. I hope all of the above shows up with proper spacing and no serious wrap around problems.
On Sep 17, 12:28&#2013266080;am, kevinjmc...@netscape.net wrote:
> On Sep 16, 4:06&#2013266080;pm, "dpjones0" <dpjon...@comcast.net> wrote: > > > > > > > Hi, > > > I have been using FFTW to do some image processing, including cross > > correlations, etc. &#2013266080;I have found it to be very useful. &#2013266080;I've used the r2c/ > > c2r combination as well as forward and reverse fftw_plan_dft_1d, based on > > examples I have found online. &#2013266080;I have been reading this discussion group > > that has helped me figure a lot of things out. &#2013266080;However, in my latest > > application, I still seem to be stuck on one issue. &#2013266080; > > > I really need to know the exact frequency values that each index in the > > FFT resultant complex represents. &#2013266080;At this point, speed is not that much of > > an issue, as much as understanding what I'm doing. &#2013266080;So for example, I set > > the input data to an NxN (e.g. 64x64) size image and use an out-of-place > > fftw_plan_dft_1d setting the input imaginary data to 0 and the real data to > > the pixel values. &#2013266080;I then end up with an NxN size FFTW complex, where the > > DC value is at index 0,0. &#2013266080;Does this also imply that in the unshifted FFTW > > frequency domain data, the first entire column has a horizontal frequency > > of 0 and the first entire row has vertical frequency of 0???. &#2013266080;I presumed > > this was the case, since the index 0,0 corresponded to 0 vertical freq and > > 0 horizontal freq, but at this point it's not clear to me. &#2013266080;Does this mean > > that for instance, along row=0, the horizontal frequency would go from: > > > (row,col) &#2013266080; &#2013266080; &#2013266080; horiz. freq. &#2013266080; &#2013266080; &#2013266080; vert. freq. > > (0,0) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 0 &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0 &#2013266080; &#2013266080;\ > > (0,1) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 1*1/N &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0 &#2013266080; &#2013266080; |-> 32 indices f=0 to > > .5-1/N > > (0,N/2 - 1) &#2013266080; &#2013266080; &#2013266080; (N/2-1)*1/N &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0 &#2013266080; &#2013266080;/ > > > (0,N/2) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; - (N/2)*1/N &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;0\ > > (0,N-1) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; - (1/(N/2) &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; 0/ &#2013266080;--> 32 indices f= -.5 to > > -1/N > > > The reason I ask this is I got confused by trying to center the DC > > component. &#2013266080;To look at power spectrums, etc with DC centered, in this > > discussion group and elsewhere online, it is suggested to move the DC > > component to N/2,N/2, and swap the quadrants (e.g. 1 and 4, 2 and 3). &#2013266080; > > > So, after doing that the DC component is now at row=N/2, col=N/2 (e.g. > > 32,32 for my 64x64 input image). &#2013266080;Then, the horizontal frequencies increase > > from coordinate row=32,col=32 to coordinate row=32,col=63 essentially from > > DC (0 freq) to PI/32 in 32 pixels. &#2013266080;However, the conjugate symmetric data > > for the negative frequencies now occurs from coordinates row=31,col=31 to > > row=31,col=1, with an extra value at row=31, col=0. &#2013266080;Why is the conjugate > > symmetric data offset by a row? &#2013266080;What are the frequencies for each of the > > points surrounding the center DC position row=32,col=32? > > > Alternatively, it was also posted that you could just multiply the input > > image pixels at coordinates (i,j) by (-1)^(i+j) for instance. &#2013266080;The problem > > that I have is that these two options do not present the same results. &#2013266080;If > > I do that, there is no DC value at N/2,N/2. &#2013266080;In other words, the data from > > row=32,col=32 to row=32,col=63 equals the conjugate symmetric data from > > row=32,col=31 to row=32,col=1, with an oddball number in the 0 column. > > > I think I'm missing something fundamental, b/c this seems a lot more > > confusing than it should be... > > > Thanks for any help in these matters. > > -Dan > > dpjones0 > > Consider an 8x8 example. &#2013266080;Consider just the FFT of one row. &#2013266080;You'll > have: > > 0f, 1f, 2f, 3f, 4f, -3f, -2f, -1f > > You go from low frequency to high frequency and then back to low > (considering -3f to be higher than -1f). &#2013266080;All of the rows will have > that frequency arrangement. &#2013266080;Now consider the FFT of one column, > You'll have: > > &#2013266080;0f > &#2013266080;1f > &#2013266080;2f > &#2013266080;3f > &#2013266080;4f > -3f > -2f > -1f > > All the columns will have the same frequency arrangement. &#2013266080;Now > consider the 4,4 point in the output - that's the (4f,4f) location - > your highest frequency. &#2013266080;The 0,0 point is DC - your lowest frequency. > The 3,3 point is (3f,3f), 2,2 is (2f,2f), etc. &#2013266080;But because of the > negative frequencies, the point 7,7 is actually (-1f,-1f). &#2013266080;Similarly, > point 6,7 is (-2f,-1f). Before centering the DC point, the 4,4 point > is the highest frequency, and you go from high frequencies to low > frequencies as you move towards the edges. > > When you center the DC output, you are in fact exchanging the data in > the 4 quadrants. &#2013266080;For N even, each one of those quadrants contains (N/ > 2)x(N/2) points. > > For the N = 8x8 case, with the DC output centered (or rotated) to 4,4, > the surrounding points are (row, column): > > (-1f, -3f) &#2013266080;(0f, -1f) &#2013266080;(1f, -1f) > (-1f, &#2013266080;0f) &#2013266080; &#2013266080;DC &#2013266080; &#2013266080; &#2013266080;(1f, 0f) > (-1f, &#2013266080;1f) &#2013266080;(0f, 1f) &#2013266080; (1f, 1f) > > OW! &#2013266080;My head hurts! &#2013266080;It really is confusing. For instance, before > centering, columns 0, 1 and 2 have the following frequency > arrangements: > > 0f0f &#2013266080; &#2013266080;1f0f &#2013266080; &#2013266080;2f0f > 0f1f &#2013266080; &#2013266080;1f1f &#2013266080; &#2013266080;2f1f > 0f2f &#2013266080; &#2013266080;1f2f &#2013266080; &#2013266080;2f2f > 0f3f &#2013266080; &#2013266080;1f3f &#2013266080; &#2013266080;2f3f > 0f4f &#2013266080; &#2013266080;1f4f &#2013266080; &#2013266080;2f4f > 0f-3f &#2013266080; 1f-3f &#2013266080; 2f-3f > 0f-2f &#2013266080; 1f-2f &#2013266080; 2f-2f > 0f-1f &#2013266080; 1f-1f &#2013266080; 2f-1f > > After centering, the arrangement gets a lot more confusing because of > the quadrant swapping. > > I've also seen the multiply outputs bit, and I think it actually does > work, because I saw some code that used it. &#2013266080;But I've always just used > 2 'for' loops to get it done. > > I hope all of the above shows up with proper spacing and no serious > wrap around problems.- Hide quoted text - > > - Show quoted text -
Oops. The point to the upper left of DC after centering is the 7,7 point, which is (-1f,-1f), not (-1f,-3f). And yes, you're right about the symmetry. Think in terms of an 8 point 1D FFT - there's no conjugate for the 0f and 4f points. And then consider the 2D case.
dpjones0 wrote:
(snip)

> I really need to know the exact frequency values that each index in the > FFT resultant complex represents. At this point, speed is not that much of > an issue, as much as understanding what I'm doing.
Just to be sure, this assumes that your data is, in fact, periodic in x and y with the period commensurate with the size of the transform...
> So for example, I set > the input data to an NxN (e.g. 64x64) size image and use an out-of-place > fftw_plan_dft_1d setting the input imaginary data to 0 and the real data to > the pixel values. I then end up with an NxN size FFTW complex, where the > DC value is at index 0,0.
The matrix axes will be x and y frequency, like components of a vector. You can then compute the frequency and direction of each by vector addition.
> Does this also imply that in the unshifted FFTW > frequency domain data, the first entire column has a horizontal frequency > of 0 and the first entire row has vertical frequency of 0???.
It should be that, just as with a 1D transform they go up from zero for the first half, and the negative frequencies for the second half. Note that FT is separable in rectangular coordinates. You can take the X FT's followed by the Y FT's, or the other way around, for the same result. -- glen