Path: news.mathworks.com!not-for-mail
From: "Chen Sagiv" <chensagivron@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Does Matlab do FFT correctly ?
Date: Sat, 3 May 2008 04:29:02 +0000 (UTC)
Organization: Sagiv Enterprises
Lines: 248
Message-ID: <fvgpme$jlv$1@fred.mathworks.com>
References: <fvfs3j$bnf$1@fred.mathworks.com> <fvg2pb$73p$1@fred.mathworks.com> <fvg3pa$slk$1@fred.mathworks.com> <fvg6sh$6nq$1@fred.mathworks.com>
Reply-To: "Chen Sagiv" <chensagivron@gmail.com>
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 1209788942 20159 172.30.248.37 (3 May 2008 04:29:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sat, 3 May 2008 04:29:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 805005
Xref: news.mathworks.com comp.soft-sys.matlab:466387



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