Thread Subject: Does Matlab do FFT correctly ?

Hi all,

We all know that the DFT of a square pulse is a sinc.

Let's try it:

% First we define the square pulse
x=-1:0.01:1;
f=zeros(size(x));
f(x>-0.5 & x<0.5)=1;
figure ; plot(x,f);

% Now let's view its Fourier transform:

figure ; plot(real(fftshift(fft(f))));
figure ; plot(imag(fftshift(fft(f))));

% Hey, the imaginary part is NOT zero
% a remedy:
figure ; plot(real(fftshift(fft(ifftshift(f)))));
figure ; plot(imag(fftshift(fft(ifftshift(f)))));

So, it seems that Matlab does correct FFT to signals that
are not centered about 0, but we have to ifftshift them
first.

Is that a known "feature" of Matlab ? I will be glad to get
a clarification !

Best,

Chen Sagiv

On May 2, 4:04 pm, "Chen Sagiv" <chensagiv...@gmail.com> wrote:
> We all know that the DFT of a square pulse is a sinc.

Actually, it's not. The *Fourier* transform of a square pulse is a
sinc function. The *discrete* Fourier transform of a square pulse is
something that is close to a sinc, but is not quite the same; it is
the ratio of two sines. (Just plug a square window into the DFT
formula and sum the geometric series.)

Your confusion is based on the fact that you don't know precisely what
you are computing.

> Let's try it:
>
> % First we define the square pulse
> x=-1:0.01:1;
> f=zeros(size(x));
> f(x>-0.5 & x<0.5)=1;
> figure ; plot(x,f);
>
> [...]
>
> So, it seems that Matlab does correct FFT to signals that
> are not centered about 0, but we have to ifftshift them
> first.

This has nothing to do with Matlab, and is a consequence of the
definition of the DFT. The origin of a DFT is at 1st element of the
input (i.e. the "left" end of the array), with periodic boundaries.
This is totally standard and more-or-less universal among
implementations of the DFT and FFT algorithms, and is not a Matlab
quirk.

The fftshift function shifts the origin to the center of the array,
which is where many people expect it (and which is often more
convenient for plotting).

Regards,
Steven G. Johnson

Dear Steven,

I appreciate your answer, but I am afraid that my question
is still open. When you run my code, you get an imaginary
part for the transform that is not supposed to be there !!!
I still have to check your remark regarding the DFT of a
square pulse being the ratio of two sines. Still, I do not
think that one may expect getting complex value for the
Fourier transform.

Best,

Chen

"Steven G. Johnson" <stevenj@alum.mit.edu> wrote in message
<80a6a285-8758-488c-9b98-
beed0d03d19e@d45g2000hsc.googlegroups.com>...
> On May 2, 4:04 pm, "Chen Sagiv" <chensagiv...@gmail.com>
wrote:
> > We all know that the DFT of a square pulse is a sinc.
>
> Actually, it's not. The *Fourier* transform of a square
pulse is a
> sinc function. The *discrete* Fourier transform of a
square pulse is
> something that is close to a sinc, but is not quite the
same; it is
> the ratio of two sines. (Just plug a square window into
the DFT
> formula and sum the geometric series.)
>
> Your confusion is based on the fact that you don't know
precisely what
> you are computing.
>
> > Let's try it:
> >
> > % First we define the square pulse
> > x=-1:0.01:1;
> > f=zeros(size(x));
> > f(x>-0.5 & x<0.5)=1;
> > figure ; plot(x,f);
> >
> > [...]
> >
> > So, it seems that Matlab does correct FFT to signals
that
> > are not centered about 0, but we have to ifftshift them
> > first.
>
> This has nothing to do with Matlab, and is a consequence
of the
> definition of the DFT. The origin of a DFT is at 1st
element of the
> input (i.e. the "left" end of the array), with periodic
boundaries.
> This is totally standard and more-or-less universal among
> implementations of the DFT and FFT algorithms, and is not
a Matlab
> quirk.
>
> The fftshift function shifts the origin to the center of
the array,
> which is where many people expect it (and which is often
more
> convenient for plotting).
>
> Regards,
> Steven G. Johnson

On May 2, 5:26 pm, "Chen Sagiv" <chensagiv...@gmail.com> wrote:
> I appreciate your answer, but I am afraid that my question
> is still open. When you run my code, you get an imaginary
> part for the transform that is not supposed to be there !!!
> I still have to check your remark regarding the DFT of a
> square pulse being the ratio of two sines. Still, I do not
> think that one may expect getting complex value for the
> Fourier transform.

Sorry, I thought I had been clear. Your square pulse is not centered
at the origin (the origin for the DFT is *not* at the centered at the
middle of the array). Because your pulse is shifted in the time
domain, it gets multiplied by a complex phase in the frequency domain
by the shift theorem.

Calling ifftshift on the fft input recenters your data correctly,
giving you a real result for the fft. However, even with this you
have to be careful, because there is some subtletly in how ifftshift
defines the "center" of the array depending upon whether the length is
even or odd.

Rest assured that Matlab's results are completely correct and
consistent with the standard definition of the DFT. Whenever it gives
you something different from what you expect, you should be looking
for a flaw in your understanding.

Regards,
Steven G. Johnson

"Chen Sagiv" <chensagivron@gmail.com> wrote in message
<fvfs3j$bnf$1@fred.mathworks.com>...
> Hi all,
>
> We all know that the DFT of a square pulse is a sinc.
>
> Let's try it:
>
> % First we define the square pulse
> x=-1:0.01:1;
> f=zeros(size(x));
> f(x>-0.5 & x<0.5)=1;
> figure ; plot(x,f);
>
> % Now let's view its Fourier transform:
>
> figure ; plot(real(fftshift(fft(f))));
> figure ; plot(imag(fftshift(fft(f))));
>
> % Hey, the imaginary part is NOT zero
> % a remedy:
> figure ; plot(real(fftshift(fft(ifftshift(f)))));
> figure ; plot(imag(fftshift(fft(ifftshift(f)))));
>
> So, it seems that Matlab does correct FFT to signals that
> are not centered about 0, but we have to ifftshift them
> first.
>
> Is that a known "feature" of Matlab ? I will be glad to
get
> a clarification !
>
> Best,
>
> Chen Sagiv
>

Chen,

Matlab does do the FFT and IFFT correctly. The problem
with your example is that your square pulse is not
perfectly symmetric. It is symmetric around zero but since
the total number of data points is odd, there has to be one
more point at -1 than +1 or the other way around.

Try generating data from -1 to (+1-stepsize) as follows,

---
% First we define the square pulse
x=-1:0.01:(1-0.01); % NOTE the change.
f=zeros(size(x));
f(x>-0.5 & x<0.5)=1;
figure ; plot(x,f);

% Now let's view its Fourier transform:

figure ; plot(real(fftshift(fft(f))));
figure ; plot(imag(fftshift(fft(f))));
---

Now the imaginary parts are very close to zero. The same
sort of behaviour happens with a sine wave that starts at
zero and ends at zero; you don't get the expected result
because the signal has a one point "glitch" if you repeat
it to + and - infinity.

Ken

Dear Steven,

Thanks, I get your point.

So, basically this is a matter of how to define the origin,
and I agree with you that it seems consistent with DFT
definition.

Nevertheless, I will have to argue that this point is not
clear at all. Most people refer only to the abs of the FT
anyway, and this problem will simply go un-noticed. Only
when phase is of interest, you really look at this subtle
point.

Best,

Chen

 "Steven G. Johnson" <stevenj@alum.mit.edu> wrote in
message <e4eb82ac-4ea7-49b7-8cdd-
aa427d8b064b@i76g2000hsf.googlegroups.com>...
> On May 2, 5:26 pm, "Chen Sagiv" <chensagiv...@gmail.com>
wrote:
> > I appreciate your answer, but I am afraid that my
question
> > is still open. When you run my code, you get an
imaginary
> > part for the transform that is not supposed to be
there !!!
> > I still have to check your remark regarding the DFT of a
> > square pulse being the ratio of two sines. Still, I do
not
> > think that one may expect getting complex value for the
> > Fourier transform.
>
> Sorry, I thought I had been clear. Your square pulse is
not centered
> at the origin (the origin for the DFT is *not* at the
centered at the
> middle of the array). Because your pulse is shifted in
the time
> domain, it gets multiplied by a complex phase in the
frequency domain
> by the shift theorem.
>
> Calling ifftshift on the fft input recenters your data
correctly,
> giving you a real result for the fft. However, even with
this you
> have to be careful, because there is some subtletly in
how ifftshift
> defines the "center" of the array depending upon whether
the length is
> even or odd.
>
> Rest assured that Matlab's results are completely correct
and
> consistent with the standard definition of the DFT.
Whenever it gives
> you something different from what you expect, you should
be looking
> for a flaw in your understanding.
>
> Regards,
> Steven G. Johnson

Dear Ken,

Thanks for your friendly answer.
I am wondering what are the implications of your input on
real life signals and images.

I am looking at the Fourier phase of signals and images.
Now, when I calculate the fft of say, an images, will I get
the correct phase, or do I have to do some tricks ?

I am not sure.

Best,

Chen

"Ken Garrard" <ken_garrardAT@ncsuDOT.edu> wrote in message
<fvg2pb$73p$1@fred.mathworks.com>...
> "Chen Sagiv" <chensagivron@gmail.com> wrote in message
> <fvfs3j$bnf$1@fred.mathworks.com>...
> > Hi all,
> >
> > We all know that the DFT of a square pulse is a sinc.
> >
> > Let's try it:
> >
> > % First we define the square pulse
> > x=-1:0.01:1;
> > f=zeros(size(x));
> > f(x>-0.5 & x<0.5)=1;
> > figure ; plot(x,f);
> >
> > % Now let's view its Fourier transform:
> >
> > figure ; plot(real(fftshift(fft(f))));
> > figure ; plot(imag(fftshift(fft(f))));
> >
> > % Hey, the imaginary part is NOT zero
> > % a remedy:
> > figure ; plot(real(fftshift(fft(ifftshift(f)))));
> > figure ; plot(imag(fftshift(fft(ifftshift(f)))));
> >
> > So, it seems that Matlab does correct FFT to signals
that
> > are not centered about 0, but we have to ifftshift them
> > first.
> >
> > Is that a known "feature" of Matlab ? I will be glad to
> get
> > a clarification !
> >
> > Best,
> >
> > Chen Sagiv
> >
>
> Chen,
>
> Matlab does do the FFT and IFFT correctly. The problem
> with your example is that your square pulse is not
> perfectly symmetric. It is symmetric around zero but
since
> the total number of data points is odd, there has to be
one
> more point at -1 than +1 or the other way around.
>
> Try generating data from -1 to (+1-stepsize) as follows,
>
> ---
> % First we define the square pulse
> x=-1:0.01:(1-0.01); % NOTE the change.
> f=zeros(size(x));
> f(x>-0.5 & x<0.5)=1;
> figure ; plot(x,f);
>
> % Now let's view its Fourier transform:
>
> figure ; plot(real(fftshift(fft(f))));
> figure ; plot(imag(fftshift(fft(f))));
> ---
>
> Now the imaginary parts are very close to zero. The same
> sort of behaviour happens with a sine wave that starts at
> zero and ends at zero; you don't get the expected result
> because the signal has a one point "glitch" if you repeat
> it to + and - infinity.
>
> Ken
>

Chen,

There are many little tricks that you can play with real
data that may help with your analysis. The best one I know
of it to use your intuition about your problem and the
answer that you expect to get. And there are several things
that I would recommend trying.
1) Look at your data carefully before doing the FFT.
   a) Center your data about its mean value
   b) Level the data to remove any drift
2) Try using windowing functions.
3) Check out www.dspguide.com for lots of practical
information on a wide variety of problems.
4) Sometimes it's helpful to do an autocovariance to see
peaks in terms of lag length (back to time and space)
instead of frequency.
5) Do exactly what you're doing - run simple examples for
which you know what the answer is supposed to be and
resolve any descrepancy. In this case the input data
wasn't exactly what it was supposed to be.

Since you asked about phase, one trick that I've found
useful is to apply a simple threshold filter to the result
returned by the FFT before proceeding with a phase angle
calculation or unwrapping process. This can also clear up
problems when doing an inverse transform and you expect a
real result and get a complex one.

A post of mine about this from a couple of years ago is
quoted below.

------------------------
I usually filter the real and imaginary parts of the
complex return vector after an FFT to force values below a
threshold to be exact zeros. Obviously these small values
have little effect on the magnitude of a signal, but they
can really foul up phase calculation.

For example, generating one period of a sine wave with 100
samples. (The last point is not a duplicate of the first)

 t=0:0.01:1-.01;
 s=sin(2*pi*t);

Then calculate
 f=fft(s);
 m=abs(f);
 p=unwrap(angle(f));

You would expect the only nonzero values in m to be at
indices 2 and 100 and that they would be 50.0. The latter
is the case; however none of the other values in m are
exactly zero.

  f(2) is -1.110223024625157e-015-50i
  f(end) is -1.110223024625157e-015+50i

and for example,
  f(33) is -1.070083830450843e-015 -6.598918101450678e-016i
  m(33) is 1.257193941000705e-015

This usually isn't a problem, although I guess it depends
on what you do next with the f or m vectors.

But the phase angles don't make much sense, wrapped or
unwrapped,
  eg, p(33) is 28.82692282385713

Again, I would expect a perfect sine wave to only have two
non-zero phase values, p(2) equal to -1.57079632679490 and
p(end) is equal to -p(2). The problem is that the inverse
tangent of a ratio is nonsense if one or both of the values
is "noise" or roundoff from a previous calculation. So we
get (in radians),

    p(2) is 4.71238898038469
and p(end) is 51.83627878423159

If the unwrap is left out you get "noise" between pi and -
pi. But if you apply a threshold filter to the result of
the FFT then the phase angles are as expected. Setting
threshold to an appropriate value say 1e-10, one way to do
this is,

function ff=cplxf(f,threshold)
 fr = real(f);
 fi = imag(f);
 irz = find(abs(fr)<threshold);
 iiz = find(abs(fi)<threshold);
 fr(irz) = 0;
 fi(iiz) = 0;
 ff = fr + i*fi;

Then,
 f=fft(s);
 ff=cplxf(f,1e-10);
 m=abs(ff);
 p=unwrap(angle(ff));

Now,
 m(2) and m(end) are 50
 all other values in m are zero
 p(2) is -1.57079632679490
 p(end) is 1.57079632679490
 and all other values in p are zero
----------------------------------

I hope this helps and good luck.
Regards,

Ken Garrard
North Carolina State University

"Chen Sagiv" <chensagivron@gmail.com> wrote in message
<fvg3pa$slk$1@fred.mathworks.com>...
> Dear Ken,
>
> Thanks for your friendly answer.
> I am wondering what are the implications of your input on
> real life signals and images.
>
> I am looking at the Fourier phase of signals and images.
> Now, when I calculate the fft of say, an images, will I
get
> the correct phase, or do I have to do some tricks ?
>
> I am not sure.
>
> Best,
>
> Chen
>
> "Ken Garrard" <ken_garrardAT@ncsuDOT.edu> wrote in
message
> <fvg2pb$73p$1@fred.mathworks.com>...
> > "Chen Sagiv" <chensagivron@gmail.com> wrote in message
> > <fvfs3j$bnf$1@fred.mathworks.com>...
> > > Hi all,
> > >
> > > We all know that the DFT of a square pulse is a sinc.
> > >
> > > Let's try it:
> > >
> > > % First we define the square pulse
> > > x=-1:0.01:1;
> > > f=zeros(size(x));
> > > f(x>-0.5 & x<0.5)=1;
> > > figure ; plot(x,f);
> > >
> > > % Now let's view its Fourier transform:
> > >
> > > figure ; plot(real(fftshift(fft(f))));
> > > figure ; plot(imag(fftshift(fft(f))));
> > >
> > > % Hey, the imaginary part is NOT zero
> > > % a remedy:
> > > figure ; plot(real(fftshift(fft(ifftshift(f)))));
> > > figure ; plot(imag(fftshift(fft(ifftshift(f)))));
> > >
> > > So, it seems that Matlab does correct FFT to signals
> that
> > > are not centered about 0, but we have to ifftshift
them
> > > first.
> > >
> > > Is that a known "feature" of Matlab ? I will be glad
to
> > get
> > > a clarification !
> > >
> > > Best,
> > >
> > > Chen Sagiv
> > >
> >
> > Chen,
> >
> > Matlab does do the FFT and IFFT correctly. The problem
> > with your example is that your square pulse is not
> > perfectly symmetric. It is symmetric around zero but
> since
> > the total number of data points is odd, there has to be
> one
> > more point at -1 than +1 or the other way around.
> >
> > Try generating data from -1 to (+1-stepsize) as follows,
> >
> > ---
> > % First we define the square pulse
> > x=-1:0.01:(1-0.01); % NOTE the change.
> > f=zeros(size(x));
> > f(x>-0.5 & x<0.5)=1;
> > figure ; plot(x,f);
> >
> > % Now let's view its Fourier transform:
> >
> > figure ; plot(real(fftshift(fft(f))));
> > figure ; plot(imag(fftshift(fft(f))));
> > ---
> >
> > Now the imaginary parts are very close to zero. The
same
> > sort of behaviour happens with a sine wave that starts
at
> > zero and ends at zero; you don't get the expected
result
> > because the signal has a one point "glitch" if you
repeat
> > it to + and - infinity.
> >
> > Ken
> >
>

Chen,

There are many little tricks that you can play with real
data that may help with your analysis. The best one I know
of it to use your intuition about your problem and the
answer that you expect to get. And there are several things
that I would recommend trying.
1) Look at your data carefully before doing the FFT.
   a) Center your data about its mean value
   b) Level the data to remove any drift
2) Try using windowing functions.
3) Check out www.dspguide.com for lots of practical
information on a wide variety of problems.
4) Sometimes it's helpful to do an autocovariance to see
peaks in terms of lag length (back to time and space)
instead of frequency.
5) Do exactly what you're doing - run simple examples for
which you know what the answer is supposed to be and
resolve any descrepancy. In this case the input data
wasn't exactly what it was supposed to be.

Since you asked about phase, one trick that I've found
useful is to apply a simple threshold filter to the result
returned by the FFT before proceeding with a phase angle
calculation or unwrapping process. This can also clear up
problems when doing an inverse transform and you expect a
real result and get a complex one.

A post of mine about this from a couple of years ago is
quoted below.

------------------------
I usually filter the real and imaginary parts of the
complex return vector after an FFT to force values below a
threshold to be exact zeros. Obviously these small values
have little effect on the magnitude of a signal, but they
can really foul up phase calculation.

For example, generating one period of a sine wave with 100
samples. (The last point is not a duplicate of the first)

 t=0:0.01:1-.01;
 s=sin(2*pi*t);

Then calculate
 f=fft(s);
 m=abs(f);
 p=unwrap(angle(f));

You would expect the only nonzero values in m to be at
indices 2 and 100 and that they would be 50.0. The latter
is the case; however none of the other values in m are
exactly zero.

  f(2) is -1.110223024625157e-015-50i
  f(end) is -1.110223024625157e-015+50i

and for example,
  f(33) is -1.070083830450843e-015 -6.598918101450678e-016i
  m(33) is 1.257193941000705e-015

This usually isn't a problem, although I guess it depends
on what you do next with the f or m vectors.

But the phase angles don't make much sense, wrapped or
unwrapped,
  eg, p(33) is 28.82692282385713

Again, I would expect a perfect sine wave to only have two
non-zero phase values, p(2) equal to -1.57079632679490 and
p(end) is equal to -p(2). The problem is that the inverse
tangent of a ratio is nonsense if one or both of the values
is "noise" or roundoff from a previous calculation. So we
get (in radians),

    p(2) is 4.71238898038469
and p(end) is 51.83627878423159

If the unwrap is left out you get "noise" between pi and -
pi. But if you apply a threshold filter to the result of
the FFT then the phase angles are as expected. Setting
threshold to an appropriate value say 1e-10, one way to do
this is,

function ff=cplxf(f,threshold)
 fr = real(f);
 fi = imag(f);
 irz = find(abs(fr)<threshold);
 iiz = find(abs(fi)<threshold);
 fr(irz) = 0;
 fi(iiz) = 0;
 ff = fr + i*fi;

Then,
 f=fft(s);
 ff=cplxf(f,1e-10);
 m=abs(ff);
 p=unwrap(angle(ff));

Now,
 m(2) and m(end) are 50
 all other values in m are zero
 p(2) is -1.57079632679490
 p(end) is 1.57079632679490
 and all other values in p are zero
----------------------------------

I hope this helps and good luck.
Regards,

Ken Garrard
North Carolina State University

"Chen Sagiv" <chensagivron@gmail.com> wrote in message
<fvg3pa$slk$1@fred.mathworks.com>...
> Dear Ken,
>
> Thanks for your friendly answer.
> I am wondering what are the implications of your input on
> real life signals and images.
>
> I am looking at the Fourier phase of signals and images.
> Now, when I calculate the fft of say, an images, will I
get
> the correct phase, or do I have to do some tricks ?
>
> I am not sure.
>
> Best,
>
> Chen
>
> "Ken Garrard" <ken_garrardAT@ncsuDOT.edu> wrote in
message
> <fvg2pb$73p$1@fred.mathworks.com>...
> > "Chen Sagiv" <chensagivron@gmail.com> wrote in message
> > <fvfs3j$bnf$1@fred.mathworks.com>...
> > > Hi all,
> > >
> > > We all know that the DFT of a square pulse is a sinc.
> > >
> > > Let's try it:
> > >
> > > % First we define the square pulse
> > > x=-1:0.01:1;
> > > f=zeros(size(x));
> > > f(x>-0.5 & x<0.5)=1;
> > > figure ; plot(x,f);
> > >
> > > % Now let's view its Fourier transform:
> > >
> > > figure ; plot(real(fftshift(fft(f))));
> > > figure ; plot(imag(fftshift(fft(f))));
> > >
> > > % Hey, the imaginary part is NOT zero
> > > % a remedy:
> > > figure ; plot(real(fftshift(fft(ifftshift(f)))));
> > > figure ; plot(imag(fftshift(fft(ifftshift(f)))));
> > >
> > > So, it seems that Matlab does correct FFT to signals
> that
> > > are not centered about 0, but we have to ifftshift
them
> > > first.
> > >
> > > Is that a known "feature" of Matlab ? I will be glad
to
> > get
> > > a clarification !
> > >
> > > Best,
> > >
> > > Chen Sagiv
> > >
> >
> > Chen,
> >
> > Matlab does do the FFT and IFFT correctly. The problem
> > with your example is that your square pulse is not
> > perfectly symmetric. It is symmetric around zero but
> since
> > the total number of data points is odd, there has to be
> one
> > more point at -1 than +1 or the other way around.
> >
> > Try generating data from -1 to (+1-stepsize) as follows,
> >
> > ---
> > % First we define the square pulse
> > x=-1:0.01:(1-0.01); % NOTE the change.
> > f=zeros(size(x));
> > f(x>-0.5 & x<0.5)=1;
> > figure ; plot(x,f);
> >
> > % Now let's view its Fourier transform:
> >
> > figure ; plot(real(fftshift(fft(f))));
> > figure ; plot(imag(fftshift(fft(f))));
> > ---
> >
> > Now the imaginary parts are very close to zero. The
same
> > sort of behaviour happens with a sine wave that starts
at
> > zero and ends at zero; you don't get the expected
result
> > because the signal has a one point "glitch" if you
repeat
> > it to + and - infinity.
> >
> > Ken
> >
>

On May 2, 6:15 pm, "Chen Sagiv" <chensagiv...@gmail.com> wrote:
> Thanks for your friendly answer.
> I am wondering what are the implications of your input on
> real life signals and images.
>
> I am looking at the Fourier phase of signals and images.
> Now, when I calculate the fft of say, an images, will I get
> the correct phase, or do I have to do some tricks ?

Your question is ill-defined, because it depends on what you mean by
"correct". From a certain point of view, the "correct" thing is what
Matlab does.

I think you mean, "will I get the phase corresponding to the origin at
the center of the image" and the answer is no, of course: you will get
the phase corresponding to the origin at the corner, corresponding to
the DFT definition, unless you use ifftshift first.

In real life, in any real application, it is easy to deal consistently
with the origin corresponding to the DFT definition, and you are never
in practice forced to use ifftshift (unless you want to). But, as
with all thing in numerical calculation, whenever you use a black box
like fft(...), you don't need to know *how* it computes what it does
but you had better understaht *what* it is computing. To do any kind
of discrete signal processing with FFTs, you have to understand the
DFT and how it differs from the continuous Fourier transform.

Regards,
Steven G. Johnson

On May 2, 6:09 pm, "Chen Sagiv" <chensagiv...@gmail.com> wrote:
> Nevertheless, I will have to argue that this point is not
> clear at all. Most people refer only to the abs of the FT
> anyway, and this problem will simply go un-noticed. Only
> when phase is of interest, you really look at this subtle
> point.

Oh, absolutely. Fourier transforms are subtle to begin with, and then
when you add discretization etc. on top of that it is a tricky
subject. Lots of students find it difficult, and digital signal
processing and DFTs are usually a relatively advanced subject in an
university curriculum. And it's easy to find people online who are
confused about this topic.

However, it's not Matlab's fault, and I don't think there's much that
they can do to eliminate the subtleties of this subject.

(If you think they should shift the origin of the FFT to the center of
the array for you, then think again. First, doing so would make them
incompatible with the universal convention of textbooks, papers, etc.
on this subject. Second, it's not a trivial matter to "center" the
origin, since the location of the "center" pixel is ambiguous for even
N; this makes the ifftshift function more tricky to use properly than
it may seem. Third, it makes the DFT formula more complicated-
looking. Fourth, the location of the origin is really the least of
the subtleties of the DFT.)

Regards,
Steven G. Johnson

Dear Ken,

Thanks for your detailed and clear answer.
I will certainly use your input.

Best,

Chen Sagiv

"Ken Garrard" <ken_garrardAT@ncsuDOT.edu> wrote in message
<fvg6sh$6nq$1@fred.mathworks.com>...
> Chen,
>
> There are many little tricks that you can play with real
> data that may help with your analysis. The best one I
know
> of it to use your intuition about your problem and the
> answer that you expect to get. And there are several
things
> that I would recommend trying.
> 1) Look at your data carefully before doing the FFT.
> a) Center your data about its mean value
> b) Level the data to remove any drift
> 2) Try using windowing functions.
> 3) Check out www.dspguide.com for lots of practical
> information on a wide variety of problems.
> 4) Sometimes it's helpful to do an autocovariance to see
> peaks in terms of lag length (back to time and space)
> instead of frequency.
> 5) Do exactly what you're doing - run simple examples for
> which you know what the answer is supposed to be and
> resolve any descrepancy. In this case the input data
> wasn't exactly what it was supposed to be.
>
> Since you asked about phase, one trick that I've found
> useful is to apply a simple threshold filter to the
result
> returned by the FFT before proceeding with a phase angle
> calculation or unwrapping process. This can also clear
up
> problems when doing an inverse transform and you expect a
> real result and get a complex one.
>
> A post of mine about this from a couple of years ago is
> quoted below.
>
> ------------------------
> I usually filter the real and imaginary parts of the
> complex return vector after an FFT to force values below
a
> threshold to be exact zeros. Obviously these small values
> have little effect on the magnitude of a signal, but they
> can really foul up phase calculation.
>
> For example, generating one period of a sine wave with
100
> samples. (The last point is not a duplicate of the
first)
>
> t=0:0.01:1-.01;
> s=sin(2*pi*t);
>
> Then calculate
> f=fft(s);
> m=abs(f);
> p=unwrap(angle(f));
>
> You would expect the only nonzero values in m to be at
> indices 2 and 100 and that they would be 50.0. The latter
> is the case; however none of the other values in m are
> exactly zero.
>
> f(2) is -1.110223024625157e-015-50i
> f(end) is -1.110223024625157e-015+50i
>
> and for example,
> f(33) is -1.070083830450843e-015 -6.598918101450678e-
016i
> m(33) is 1.257193941000705e-015
>
> This usually isn't a problem, although I guess it depends
> on what you do next with the f or m vectors.
>
> But the phase angles don't make much sense, wrapped or
> unwrapped,
> eg, p(33) is 28.82692282385713
>
> Again, I would expect a perfect sine wave to only have
two
> non-zero phase values, p(2) equal to -1.57079632679490
and
> p(end) is equal to -p(2). The problem is that the inverse
> tangent of a ratio is nonsense if one or both of the
values
> is "noise" or roundoff from a previous calculation. So we
> get (in radians),
>
> p(2) is 4.71238898038469
> and p(end) is 51.83627878423159
>
> If the unwrap is left out you get "noise" between pi and -
> pi. But if you apply a threshold filter to the result of
> the FFT then the phase angles are as expected. Setting
> threshold to an appropriate value say 1e-10, one way to
do
> this is,
>
> function ff=cplxf(f,threshold)
> fr = real(f);
> fi = imag(f);
> irz = find(abs(fr)<threshold);
> iiz = find(abs(fi)<threshold);
> fr(irz) = 0;
> fi(iiz) = 0;
> ff = fr + i*fi;
>
> Then,
> f=fft(s);
> ff=cplxf(f,1e-10);
> m=abs(ff);
> p=unwrap(angle(ff));
>
> Now,
> m(2) and m(end) are 50
> all other values in m are zero
> p(2) is -1.57079632679490
> p(end) is 1.57079632679490
> and all other values in p are zero
> ----------------------------------
>
> I hope this helps and good luck.
> Regards,
>
> Ken Garrard
> North Carolina State University
>
> "Chen Sagiv" <chensagivron@gmail.com> wrote in message
> <fvg3pa$slk$1@fred.mathworks.com>...
> > Dear Ken,
> >
> > Thanks for your friendly answer.
> > I am wondering what are the implications of your input
on
> > real life signals and images.
> >
> > I am looking at the Fourier phase of signals and
images.
> > Now, when I calculate the fft of say, an images, will I
> get
> > the correct phase, or do I have to do some tricks ?
> >
> > I am not sure.
> >
> > Best,
> >
> > Chen
> >
> > "Ken Garrard" <ken_garrardAT@ncsuDOT.edu> wrote in
> message
> > <fvg2pb$73p$1@fred.mathworks.com>...
> > > "Chen Sagiv" <chensagivron@gmail.com> wrote in
message
> > > <fvfs3j$bnf$1@fred.mathworks.com>...
> > > > Hi all,
> > > >
> > > > We all know that the DFT of a square pulse is a
sinc.
> > > >
> > > > Let's try it:
> > > >
> > > > % First we define the square pulse
> > > > x=-1:0.01:1;
> > > > f=zeros(size(x));
> > > > f(x>-0.5 & x<0.5)=1;
> > > > figure ; plot(x,f);
> > > >
> > > > % Now let's view its Fourier transform:
> > > >
> > > > figure ; plot(real(fftshift(fft(f))));
> > > > figure ; plot(imag(fftshift(fft(f))));
> > > >
> > > > % Hey, the imaginary part is NOT zero
> > > > % a remedy:
> > > > figure ; plot(real(fftshift(fft(ifftshift(f)))));
> > > > figure ; plot(imag(fftshift(fft(ifftshift(f)))));
> > > >
> > > > So, it seems that Matlab does correct FFT to
signals
> > that
> > > > are not centered about 0, but we have to ifftshift
> them
> > > > first.
> > > >
> > > > Is that a known "feature" of Matlab ? I will be
glad
> to
> > > get
> > > > a clarification !
> > > >
> > > > Best,
> > > >
> > > > Chen Sagiv
> > > >
> > >
> > > Chen,
> > >
> > > Matlab does do the FFT and IFFT correctly. The
problem
> > > with your example is that your square pulse is not
> > > perfectly symmetric. It is symmetric around zero but
> > since
> > > the total number of data points is odd, there has to
be
> > one
> > > more point at -1 than +1 or the other way around.
> > >
> > > Try generating data from -1 to (+1-stepsize) as
follows,
> > >
> > > ---
> > > % First we define the square pulse
> > > x=-1:0.01:(1-0.01); % NOTE the change.
> > > f=zeros(size(x));
> > > f(x>-0.5 & x<0.5)=1;
> > > figure ; plot(x,f);
> > >
> > > % Now let's view its Fourier transform:
> > >
> > > figure ; plot(real(fftshift(fft(f))));
> > > figure ; plot(imag(fftshift(fft(f))));
> > > ---
> > >
> > > Now the imaginary parts are very close to zero. The
> same
> > > sort of behaviour happens with a sine wave that
starts
> at
> > > zero and ends at zero; you don't get the expected
> result
> > > because the signal has a one point "glitch" if you
> repeat
> > > it to + and - infinity.
> > >
> > > Ken
> > >
> >
>

Dear Steven,

Thanks for your answer.

Best,

Chen

"Steven G. Johnson" <stevenj@alum.mit.edu> wrote in message
<da5827de-b6d8-4f66-9943-
86aa762b2175@8g2000hse.googlegroups.com>...
> On May 2, 6:15 pm, "Chen Sagiv" <chensagiv...@gmail.com>
wrote:
> > Thanks for your friendly answer.
> > I am wondering what are the implications of your input
on
> > real life signals and images.
> >
> > I am looking at the Fourier phase of signals and images.
> > Now, when I calculate the fft of say, an images, will I
get
> > the correct phase, or do I have to do some tricks ?
>
> Your question is ill-defined, because it depends on what
you mean by
> "correct". From a certain point of view, the "correct"
thing is what
> Matlab does.
>
> I think you mean, "will I get the phase corresponding to
the origin at
> the center of the image" and the answer is no, of course:
you will get
> the phase corresponding to the origin at the corner,
corresponding to
> the DFT definition, unless you use ifftshift first.
>
> In real life, in any real application, it is easy to deal
consistently
> with the origin corresponding to the DFT definition, and
you are never
> in practice forced to use ifftshift (unless you want
to). But, as
> with all thing in numerical calculation, whenever you use
a black box
> like fft(...), you don't need to know *how* it computes
what it does
> but you had better understaht *what* it is computing. To
do any kind
> of discrete signal processing with FFTs, you have to
understand the
> DFT and how it differs from the continuous Fourier
transform.
>
> Regards,
> Steven G. Johnson

Hi again,

I am not here to blame Matlab for anything, but to
understand better what Matlab is doing.

I think that a short white paper would do the work. If you
know one, I will be happy to get a pointer. If not, I will
try to find the time to write something up.

In any case, I am very happy with the fruitful discussion
that took place here, and I certainly learnt quite a lot.

Best and thanks for your efforts.

Chen

"Steven G. Johnson" <stevenj@alum.mit.edu> wrote in message
<3d638b38-6307-4a73-b305-
b1eda3ebbbbd@34g2000hsh.googlegroups.com>...
> On May 2, 6:09 pm, "Chen Sagiv" <chensagiv...@gmail.com>
wrote:
> > Nevertheless, I will have to argue that this point is
not
> > clear at all. Most people refer only to the abs of the
FT
> > anyway, and this problem will simply go un-noticed. Only
> > when phase is of interest, you really look at this
subtle
> > point.
>
> Oh, absolutely. Fourier transforms are subtle to begin
with, and then
> when you add discretization etc. on top of that it is a
tricky
> subject. Lots of students find it difficult, and digital
signal
> processing and DFTs are usually a relatively advanced
subject in an
> university curriculum. And it's easy to find people
online who are
> confused about this topic.
>
> However, it's not Matlab's fault, and I don't think
there's much that
> they can do to eliminate the subtleties of this subject.
>
> (If you think they should shift the origin of the FFT to
the center of
> the array for you, then think again. First, doing so
would make them
> incompatible with the universal convention of textbooks,
papers, etc.
> on this subject. Second, it's not a trivial matter
to "center" the
> origin, since the location of the "center" pixel is
ambiguous for even
> N; this makes the ifftshift function more tricky to use
properly than
> it may seem. Third, it makes the DFT formula more
complicated-
> looking. Fourth, the location of the origin is really
the least of
> the subtleties of the DFT.)
>
> Regards,
> Steven G. Johnson