Reply by Steven G. Johnson●August 5, 200820080805
On Aug 5, 12:05 pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> The documentation of a technical library should not attempt
> to be a tutorial. However, it might be useful for users if
> the documentation provides pointers to good tutorials.
It actually does provide a pointer to at least one good textbook on
applications of FFTs, and the links page lists many more resources
online.
Reply by Rune Allnor●August 5, 200820080805
On 5 Aug, 03:59, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> However, as the introduction to the manual explains, the manual makes
> no attempt to be a tutorial on the meaning and properties of the DFT.
> There are plenty of textbooks on that. The manual assumes that you
> know what DFT (or DCT, or ...) you want to compute and why, and its
> only function is to tell you how to use FFTW to perform that
> computation.
The documentation of a technical library should not attempt
to be a tutorial. However, it might be useful for users if
the documentation provides pointers to good tutorials.
Rune
Reply by Steven G. Johnson●August 4, 200820080804
On Aug 3, 9:17 am, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> I should also probably mention that I'm implementing that using FFTW 3
> (no Matlab to toy around with), and that its documentation is pretty
> poor on the very topic of my problem. It seemingly doesn't tell what
> the output matches to depending on the transform. For example I
> noticed that the output is quite different between a 2D Discrete
> Hartley Transform and a 2D DFT (and also a 2D DCT, but I know why that
> one is different, I at least understand that, but the documentation
> won't say anything about it),
Note that the section "what FFTW really computes" of the manual gives
the exact formulas for the transforms that it performs.
However, as the introduction to the manual explains, the manual makes
no attempt to be a tutorial on the meaning and properties of the DFT.
There are plenty of textbooks on that. The manual assumes that you
know what DFT (or DCT, or ...) you want to compute and why, and its
only function is to tell you how to use FFTW to perform that
computation.
Regards,
Steven G. Johnson
Reply by JimAtQuarktet●August 4, 200820080804
On Aug 3, 9:17 am, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> Here's my problem. I have animageto be processed (well, deconvolved
> pretty much), and a starshaped impulse response I want to deconvolve
> it with. The impulse response is symmetric in both the vertical and
> horizontal axis and pretty much centred, but not very precisely. Of
> course both theimageand the impulse response are appropriately zero
> padded as to avoid any circular convolution.
>
On the PSF topic, you (and others) may want to look at Tria at
http://www.quarktet.com/Tria.html. It has builtin centering and
cleaning functions for PSF. If your PSF is not centered before
deconvolution, your processed image will shift by the difference
between the center of the image and the center of the PSF. This may
or may not be important. Tria also uses a noniterative pseudo
inverse filter (PIF) (sometimes called the Wiener filter) that,a s you
can see from the image gallery does not produce edge effects when
applied correctly. The ImageJ plugins, as far as I can tell, are
iterative processes that some people feel are more effective than the
PIF. However, there are a few papers that demonstrate that this is
not the case. Tria is free to try, and costs more than an image
processing book, but much less than MatLab.
Best Regards,
Jim C
Reply by glen herrmannsfeldt●August 4, 200820080804
Michel Rouzic wrote:
(snip)
> Instead, it seems that 2D DFTs are achieved by performing 1D DFTs on
> each row, and then performing 1D DFTs on each column. If I understand
> correctly, the result is a complex array representing from Nyquist to
> +Nyquist on each axis, and is symmetrical on one axis due to the fact
> that the first series of 1D DFTs is performed on real data, as the
> rest is on complex data and therefore has nonredundant negative
> frequency components. My problem with that is that I don't understand
> why we should "privilege" an axis by processing it first, and
> therefore obtain a different result that if we have processed the
> other axis first.
I am pretty sure that it doesn't depend on the order of the
transforms. The Fourier transform is a linear operator.
One way to look at it is through the normal modes of a
square drum head. (Well, that would by 2D DST, but the
symmetry is the same.)
In the 2D transform array the upper left corner is the DC component,
or, for a drum head, the mode with no nodal lines. As you go across
the matrix in one direction you increase the number of vertical node
lines. In the other direction horizontal node lines.
Do the transform by hand for a 2x2 or 3x3 matrix. You will
then see the symmetry.
 glen
Reply by Chris Bore●August 4, 200820080804
You might look up ImageJ, which is public domain image processing
software and has facilities to directly manipulate 2D FFT data. It
also, by the way, has several implementations of deconvolution (as
plugins).
Chris
==================
Chris Bore
BORES Signal Processing
www.bores.com
On Aug 3, 2:17�pm, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> Here's my problem. I have an image to be processed (well, deconvolved
> pretty much), and a starshaped impulse response I want to deconvolve
> it with. The impulse response is symmetric in both the vertical and
> horizontal axis and pretty much centred, but not very precisely. Of
> course both the image and the impulse response are appropriately zero
> padded as to avoid any circular convolution.
>
> My idea is to transform the IR into a frequencydomain (magnitude
> only) representation using something like a 2D FFT, hard limit the
> frequencydomain magnitude bins under a threshold value to make sure
> that the weaker frequency components don't get amplified more than
> they should, and then perform the deconvolution.
>
> Well at this point I guess I could simply do a f(x) = 1/x on all the
> magnitude bins, zero the phase (as to make the IR centered around 0,
> 0), IFFT and convolve in the spacedomain, or I could FFT the image to
> be deconvolved and multiply (or rather divide) with the frequency
> domain magnitude of the hardlimited impulse response. The problem is
> that I don't understand how to get the magnitude of the impulse
> response FFT (or similar).
>
> There's something I don't understand about 2D DFTs. Just like the DFT
> of a 1D real signal is represented by an array of complex bins
> representing from DC to the Nyquist frequency (omitting the redundant
> negative frequencies) I would expect a 2D DFT of a (real) image to be
> a complex array ranging from each axis' DC to each axis' Nyquist
> frequency (although I realise this is surely not possible since it
> would take twice less space to represent a real image under that form)
>
> Instead, it seems that 2D DFTs are achieved by performing 1D DFTs on
> each row, and then performing 1D DFTs on each column. If I understand
> correctly, the result is a complex array representing from Nyquist to
> +Nyquist on each axis, and is symmetrical on one axis due to the fact
> that the first series of 1D DFTs is performed on real data, as the
> rest is on complex data and therefore has nonredundant negative
> frequency components. My problem with that is that I don't understand
> why we should "privilege" an axis by processing it first, and
> therefore obtain a different result that if we have processed the
> other axis first.
>
> Also, I don't understand what to do with that. If I'm not mistaken, to
> obtain the magnitude of my IR in the frequency domain, I just need to
> do a sqrt(real� + imaginary�) on each of my complex bins, correct? So
> then I'm left with.. a magnitude that is nonsymmetrical around DC on
> only one axis? What do I do with that, I mean, why should it have
> negative frequency components that are not the same as the positive
> ones, and most importantly, why on only one (arbitrarily chose) axis??
>
> I should also probably mention that I'm implementing that using FFTW 3
> (no Matlab to toy around with), and that its documentation is pretty
> poor on the very topic of my problem. It seemingly doesn't tell what
> the output matches to depending on the transform. For example I
> noticed that the output is quite different between a 2D Discrete
> Hartley Transform and a 2D DFT (and also a 2D DCT, but I know why that
> one is different, I at least understand that, but the documentation
> won't say anything about it), why? Disclaimer : I don't really know
> what the DHT does, but I find it striking that its output (as a real
> image, so I take it it's in their "half complex" format that goes DC
> Re1 Re2 Re3 �Nyquist Im3 Im2 Im1) is made of only two parts (which I
> guess would come as one if it were to be put as a complex array as it
> may (or may not? the documentation doesn't seem to discuss this) be in
> "half complex" format).
>
> Thanks in advance
Reply by Gordon Sande●August 3, 200820080803
On 20080803 14:29:27 0300, Michel Rouzic <Michel0528@yahoo.fr> said:
>
>
> Rune Allnor wrote:
>> If we start with the computations, then you are right in that
>> the concept is to first perform 1D FFTs in one direction and
>> then 1D FFTs in the other direction (although it may be possible
>> to implement the 2d FFT in other ways; I don't know).
>>
>> The end result is, however, the same regardless of the order the
>> FFTs are computed.
>
> Ha, I don't understand that part. How can it be the same thing
> considered that, if we start from a realonly 2D image, and that
> therefore the result should be symmetrical for the first processed
> axis?
No! Only along y=0. A special case only of the inversion symmetry
about the origin. It is also symmetric along y for x=0. Again another
special case.
To lower your confusion you should drop the "optimization" of using
the real values only. It has confused you no end. Work it through
without the "optimization" and only after it is "working" try to
follow the effects of real values. I use "working" as deconvolution
is a noise enhancing process. You will almost surely end up needing
some form of regualrization.
>> As for the structure of the spectrum, it will become a little bit
>> different than you expect. Assume, for the sake of argument, that
>> you compute the FFTs in the x direction first:
>>
>> f(x,y) [FFT_x]> f(u,y)
>>
>> Here f(x,y) is the realvalued image, FFT_x is the FFT along
>> the x dimension, and u is the 'spatial x frequency'.
>>
>> Since the image f(x,y) is realvalued, each line f(u) in the
>> image is complexvalued and conjugate symmetric. No surprises.
>>
>> Now, once we compute the FFT in the other dimension,
>>
>> f(u,y) [FFT_y]> f(u,v)
>>
>> the input to the ydirection FFT, the columns of f(u,y),
>> are complexvalued, not realvalued. This means that the
>> columns f(v) in the end spectrum are no longer symmetric.
>> Remember, Nyquist's theorem
>>
>> f < Fs/2
>>
>> only applies to realvalued signals because of the
>> conjugate symmetry of Euler's equations,
>>
>> cos(x) = 1/2 (exp(jx)+exp(jx))
>> sin(x) = 1/j2(exp(jx)exp(jx))
>>
>> (are the signs correct here? i don't remember...)
>>
>> Since complexvalued signals don't have the conjugate
>> symmetric spectrum, the complex valued form of Nyquist's
>> theorem becomes
>>
>> f < Fs.
>>
>> The net effect is that the 2D spectrum is cylindrical
>> conjugate symmetric around origo (view with fixedwidth
>> font):
>> ^ v
>> 
>> +++ Fy/2
>>   
>>   
>>   x (u,v) 
>> 0 +++>
>>  x   u
>>  (u,v)  
>>   
>> +++ Fy/2
>> Fx/2 0 Fx/2
>>
>> The conjugate mirror of a sinusoidal at (u,v) will
>> occur at (u.v) in a 2D spectrum centered at (0,0).
>>
>>> Also, I don't understand what to do with that. If I'm not mistaken, to
>>> obtain the magnitude of my IR in the frequency domain, I just need to
>>> do a sqrt(real� + imaginary�) on each of my complex bin
> s, correct? So
>>> then I'm left with.. a magnitude that is nonsymmetrical around DC on
>>> only one axis? What do I do with that, I mean, why should it have
>>> negative frequency components that are not the same as the positive
>>> ones, and most importantly, why on only one (arbitrarily chose) axis??
>>
>> These things are explained in some detail in the Gonzales and Woods
>> book.
>> The trick is to use circular symmetry around origo.
>
> Ha. thanks for the explanation, that clears up a few things about the
> theoretical background. But obviously I'm not going to order a $120
> book to solve this particular problem (although I'll surely consider
> the book again if it appears to cover the image processing problems I
> have yet to cover) so would you mind elaborating a bit on that please?
>
> Here's what I figure should be the way to do it, based on your
> explanations of 2D FFTs : I should do my hard limiting on all
> magnitude bins alike regardless of which frequencies they represent,
> then divide the image's FFT's complex bins with the aforementioned
> magnitude. Considered both arrays would be the same size to begin with
> (since they would appropriately padded) their spatial frequencies will
> map so I don't even really have to wonder what I'm looking at and
> simply process bins all the same. Does it sound correct?
Reply by Rune Allnor●August 3, 200820080803
On 3 Aug, 19:29, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> Rune Allnor wrote:
> > If we start with the computations, then you are right in that
> > the concept is to first perform 1D FFTs in one direction and
> > then 1D FFTs in the other direction (although it may be possible
> > to implement the 2d FFT in other ways; I don't know).
>
> > The end result is, however, the same regardless of the order the
> > FFTs are computed.
>
> Ha, I don't understand that part. How can it be the same thing
> considered that, if we start from a realonly 2D image, and that
> therefore the result should be symmetrical for the first processed
> axis?
The 2D DFT rotational symmetric around origo, not around any line
in (u,v) domain.
> > As for the structure of the spectrum, it will become a little bit
> > different than you expect. Assume, for the sake of argument, that
> > you compute the FFTs in the x direction first:
>
> > f(x,y) [FFT_x]> f(u,y)
>
> > Here f(x,y) is the realvalued image, FFT_x is the FFT along
> > the x dimension, and u is the 'spatial x frequency'.
>
> > Since the image f(x,y) is realvalued, each line f(u) in the
> > image is complexvalued and conjugate symmetric. No surprises.
>
> > Now, once we compute the FFT in the other dimension,
>
> > f(u,y) [FFT_y]> f(u,v)
>
> > the input to the ydirection FFT, the columns of f(u,y),
> > are complexvalued, not realvalued. This means that the
> > columns f(v) in the end spectrum are no longer symmetric.
> > Remember, Nyquist's theorem
>
> > f < Fs/2
>
> > only applies to realvalued signals because of the
> > conjugate symmetry of Euler's equations,
>
> > cos(x) = 1/2 (exp(jx)+exp(jx))
> > sin(x) = 1/j2(exp(jx)exp(jx))
>
> > (are the signs correct here? i don't remember...)
>
> > Since complexvalued signals don't have the conjugate
> > symmetric spectrum, the complex valued form of Nyquist's
> > theorem becomes
>
> > f < Fs.
>
> > The net effect is that the 2D spectrum is cylindrical
> > conjugate symmetric around origo (view with fixedwidth
> > font):
> > � � � � � � � � � � �^ v
> > � � � � � � � � � � �
> > � � +++ Fy/2
> > � �  � � � � � � � � � � � � � � � �
> > � �  � � � � � � � � � � � � � � � �
> > � �  � � � � � � � � � �x (u,v) � � 
> > �0 +++>
> > � �  � � � � �x � �  � � � � � � � � �u
> > � �  � � �(u,v) �  � � � � � � � �
> > � �  � � � � � � � � � � � � � � � �
> > � � +++ Fy/2
> > � Fx/2 � � � � � � �0 � � � � � � � Fx/2
>
> > The conjugate mirror of a sinusoidal at (u,v) will
> > occur at (u.v) in a 2D spectrum centered at (0,0).
>
> > > Also, I don't understand what to do with that. If I'm not mistaken, to
> > > obtain the magnitude of my IR in the frequency domain, I just need to
> > > do a sqrt(real + imaginary ) on each of my complex bins, correct? So
> > > then I'm left with.. a magnitude that is nonsymmetrical around DC on
> > > only one axis? What do I do with that, I mean, why should it have
> > > negative frequency components that are not the same as the positive
> > > ones, and most importantly, why on only one (arbitrarily chose) axis??
>
> > These things are explained in some detail in the Gonzales and Woods
> > book.
> > The trick is to use circular symmetry around origo.
>
> Ha. thanks for the explanation, that clears up a few things about the
> theoretical background. But obviously I'm not going to order a $120
> book to solve this particular problem (although I'll surely consider
> the book again if it appears to cover the image processing problems I
> have yet to cover) so would you mind elaborating a bit on that please?
I can't. It's one of those things that has become second nature
to me; it's so obvious in my mind (I know, it's not to others,
but I spent 6 years working with 2D spectra every day), so I just
don't know what to explain.
> Here's what I figure should be the way to do it, based on your
> explanations of 2D FFTs : I should do my hard limiting on all
> magnitude bins alike regardless of which frequencies they represent,
> then divide the image's FFT's complex bins with the aforementioned
> magnitude. Considered both arrays would be the same size to begin with
> (since they would appropriately padded) their spatial frequencies will
> map so I don't even really have to wonder what I'm looking at and
> simply process bins all the same. Does it sound correct?
The algorithm seems more or less correct. It's just that the $120
book shows a compariion between that approach and a few others
where this gives the by far worst result.
If possoble, pop into your local university bookstore to browse
a copy. Depending on what you want to do with the processing, it
might just save you a lot of time and grief.
Rune
Reply by Michel Rouzic●August 3, 200820080803
Rune Allnor wrote:
> If we start with the computations, then you are right in that
> the concept is to first perform 1D FFTs in one direction and
> then 1D FFTs in the other direction (although it may be possible
> to implement the 2d FFT in other ways; I don't know).
>
> The end result is, however, the same regardless of the order the
> FFTs are computed.
Ha, I don't understand that part. How can it be the same thing
considered that, if we start from a realonly 2D image, and that
therefore the result should be symmetrical for the first processed
axis?
> As for the structure of the spectrum, it will become a little bit
> different than you expect. Assume, for the sake of argument, that
> you compute the FFTs in the x direction first:
>
> f(x,y) [FFT_x]> f(u,y)
>
> Here f(x,y) is the realvalued image, FFT_x is the FFT along
> the x dimension, and u is the 'spatial x frequency'.
>
> Since the image f(x,y) is realvalued, each line f(u) in the
> image is complexvalued and conjugate symmetric. No surprises.
>
> Now, once we compute the FFT in the other dimension,
>
> f(u,y) [FFT_y]> f(u,v)
>
> the input to the ydirection FFT, the columns of f(u,y),
> are complexvalued, not realvalued. This means that the
> columns f(v) in the end spectrum are no longer symmetric.
> Remember, Nyquist's theorem
>
> f < Fs/2
>
> only applies to realvalued signals because of the
> conjugate symmetry of Euler's equations,
>
> cos(x) = 1/2 (exp(jx)+exp(jx))
> sin(x) = 1/j2(exp(jx)exp(jx))
>
> (are the signs correct here? i don't remember...)
>
> Since complexvalued signals don't have the conjugate
> symmetric spectrum, the complex valued form of Nyquist's
> theorem becomes
>
> f < Fs.
>
> The net effect is that the 2D spectrum is cylindrical
> conjugate symmetric around origo (view with fixedwidth
> font):
> ^ v
> 
> +++ Fy/2
>   
>   
>   x (u,v) 
> 0 +++>
>  x   u
>  (u,v)  
>   
> +++ Fy/2
> Fx/2 0 Fx/2
>
> The conjugate mirror of a sinusoidal at (u,v) will
> occur at (u.v) in a 2D spectrum centered at (0,0).
>
> > Also, I don't understand what to do with that. If I'm not mistaken, to
> > obtain the magnitude of my IR in the frequency domain, I just need to
> > do a sqrt(real� + imaginary�) on each of my complex bins, correct? So
> > then I'm left with.. a magnitude that is nonsymmetrical around DC on
> > only one axis? What do I do with that, I mean, why should it have
> > negative frequency components that are not the same as the positive
> > ones, and most importantly, why on only one (arbitrarily chose) axis??
>
> These things are explained in some detail in the Gonzales and Woods
> book.
> The trick is to use circular symmetry around origo.
Ha. thanks for the explanation, that clears up a few things about the
theoretical background. But obviously I'm not going to order a $120
book to solve this particular problem (although I'll surely consider
the book again if it appears to cover the image processing problems I
have yet to cover) so would you mind elaborating a bit on that please?
Here's what I figure should be the way to do it, based on your
explanations of 2D FFTs : I should do my hard limiting on all
magnitude bins alike regardless of which frequencies they represent,
then divide the image's FFT's complex bins with the aforementioned
magnitude. Considered both arrays would be the same size to begin with
(since they would appropriately padded) their spatial frequencies will
map so I don't even really have to wonder what I'm looking at and
simply process bins all the same. Does it sound correct?
Reply by Rune Allnor●August 3, 200820080803
On 3 Aug, 15:17, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> Here's my problem. I have an image to be processed (well, deconvolved
> pretty much), and a starshaped impulse response I want to deconvolve
> it with. The impulse response is symmetric in both the vertical and
> horizontal axis and pretty much centred, but not very precisely. Of
> course both the image and the impulse response are appropriately zero
> padded as to avoid any circular convolution.
>
> My idea is to transform the IR into a frequencydomain (magnitude
> only) representation using something like a 2D FFT, hard limit the
> frequencydomain magnitude bins under a threshold value to make sure
> that the weaker frequency components don't get amplified more than
> they should, and then perform the deconvolution.
Deconvolution is tricky business at the best of times. You might
want to consult the 2008 edition of the Gozales and Woods book
on image processing to see a number of approaches and some
comparisions of results.
> Well at this point I guess I could simply do a f(x) = 1/x on all the
> magnitude bins, zero the phase (as to make the IR centered around 0,
> 0), IFFT and convolve in the spacedomain, or I could FFT the image to
> be deconvolved and multiply (or rather divide) with the frequency
> domain magnitude of the hardlimited impulse response. The problem is
> that I don't understand how to get the magnitude of the impulse
> response FFT (or similar).
OK. The 2D FFT is the problem here. Let's focus on that and forget
about deconvolution for now.
> There's something I don't understand about 2D DFTs. Just like the DFT
> of a 1D real signal is represented by an array of complex bins
> representing from DC to the Nyquist frequency (omitting the redundant
> negative frequencies) I would expect a 2D DFT of a (real) image to be
> a complex array ranging from each axis' DC to each axis' Nyquist
> frequency (although I realise this is surely not possible since it
> would take twice less space to represent a real image under that form)
>
> Instead, it seems that 2D DFTs are achieved by performing 1D DFTs on
> each row, and then performing 1D DFTs on each column. If I understand
> correctly, the result is a complex array representing from Nyquist to
> +Nyquist on each axis, and is symmetrical on one axis due to the fact
> that the first series of 1D DFTs is performed on real data, as the
> rest is on complex data and therefore has nonredundant negative
> frequency components. My problem with that is that I don't understand
> why we should "privilege" an axis by processing it first, and
> therefore obtain a different result that if we have processed the
> other axis first.
If we start with the computations, then you are right in that
the concept is to first perform 1D FFTs in one direction and
then 1D FFTs in the other direction (although it may be possible
to implement the 2d FFT in other ways; I don't know).
The end result is, however, the same regardless of the order the
FFTs are computed.
As for the structure of the spectrum, it will become a little bit
different than you expect. Assume, for the sake of argument, that
you compute the FFTs in the x direction first:
f(x,y) [FFT_x]> f(u,y)
Here f(x,y) is the realvalued image, FFT_x is the FFT along
the x dimension, and u is the 'spatial x frequency'.
Since the image f(x,y) is realvalued, each line f(u) in the
image is complexvalued and conjugate symmetric. No surprises.
Now, once we compute the FFT in the other dimension,
f(u,y) [FFT_y]> f(u,v)
the input to the ydirection FFT, the columns of f(u,y),
are complexvalued, not realvalued. This means that the
columns f(v) in the end spectrum are no longer symmetric.
Remember, Nyquist's theorem
f < Fs/2
only applies to realvalued signals because of the
conjugate symmetry of Euler's equations,
cos(x) = 1/2 (exp(jx)+exp(jx))
sin(x) = 1/j2(exp(jx)exp(jx))
(are the signs correct here? i don't remember...)
Since complexvalued signals don't have the conjugate
symmetric spectrum, the complex valued form of Nyquist's
theorem becomes
f < Fs.
The net effect is that the 2D spectrum is cylindrical
conjugate symmetric around origo (view with fixedwidth
font):
^ v

+++ Fy/2
  
  
  x (u,v) 
0 +++>
 x   u
 (u,v)  
  
+++ Fy/2
Fx/2 0 Fx/2
The conjugate mirror of a sinusoidal at (u,v) will
occur at (u.v) in a 2D spectrum centered at (0,0).
> Also, I don't understand what to do with that. If I'm not mistaken, to
> obtain the magnitude of my IR in the frequency domain, I just need to
> do a sqrt(real� + imaginary�) on each of my complex bins, correct? So
> then I'm left with.. a magnitude that is nonsymmetrical around DC on
> only one axis? What do I do with that, I mean, why should it have
> negative frequency components that are not the same as the positive
> ones, and most importantly, why on only one (arbitrarily chose) axis??
These things are explained in some detail in the Gonzales and Woods
book.
The trick is to use circular symmetry around origo.
Rune