Path: news.mathworks.com!not-for-mail
From: "Ken Garrard" <ken_garrardAT@ncsuDOT.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Does Matlab do FFT correctly ?
Date: Fri, 2 May 2008 23:08:02 +0000 (UTC)
Organization: North Carolina State University
Lines: 214
Message-ID: <fvg6si$6nt$1@fred.mathworks.com>
References: <fvfs3j$bnf$1@fred.mathworks.com> <fvg2pb$73p$1@fred.mathworks.com> <fvg3pa$slk$1@fred.mathworks.com>
Reply-To: "Ken Garrard" <ken_garrardAT@ncsuDOT.edu>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1209769682 6909 172.30.248.37 (2 May 2008 23:08:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 2 May 2008 23:08:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 216874
Xref: news.mathworks.com comp.soft-sys.matlab:466372


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
> > 
>