nonzero phase spectrum from FFT of even function with zero imag part
in [DSP]
|
Prev: Phase-frequency detector with synchronous output?
Next: nonzero phase spectrum from FFT of even function with zero imag part
From: Clay on 25 Jun 2010 10:21 On Jun 24, 6:15 pm, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> wrote: > >On Jun 24, 2:12=A0pm, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> > >wrote: > >... > >> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should > gi= > >ve > >> zero.... > > >> why not? > > >> thanks > >> fisico32 > > >When I try: > >>> x =3D fft([1 1 1 1 1 1]) > > >x =3D > > > 6 0 0 0 0 0 > > >>> atan2(imag(x),real(x)) > > >ans =3D > > > 0 0 0 0 0 0 > >It works. > > >Do you have an example with a problem? In fact, when you discover > >magic like this would you always include the example in your initial > >post please. > > >Dale B. Dalrymple > > thanks Dale, > here what matlab does: > > f=[0 1 2 3 3 2 1]; (even sequence). > A=fft(f); > B=imag(A); (it is all zeros). > C=angle(A) (it shows zero and pi values...) > > The function angle should just implement the inverse tangent...- Hide quoted text - > > - Show quoted text - You know that arg(0+0j) is undefined -- right? Its a form of 0/0. So when either component has some garbage and is not identically zero then you can see values like pi pop up. Sometimes programmers put in things to try to force such results to a consistant 0 value. Clay
From: Fred Marshall on 25 Jun 2010 14:56 fisico32 wrote: >> On Jun 24, 5:12 pm, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> >> wrote: >>> Hello Forum, >>> >>> the fft of a rectangular, even, symmetric sequence show a nonzero real > part >>> and a zero imaginary part. All good so far. >>> >>> However the phase spectrum is nonzero: zero at some frequencies, pi at > some >>> frequencies , -pi at some others.... >>> >>> The phase spectrum is given by atan(Imag(fft)/real(fft)), >> Incorrect. Look at the source code: Enter into the command line. >> >> type angle >> type atan2 >> >>> so it should give zero.... >>> >>> why not? >> 1. Choose your own examples and compare atan(y/x) with atan2(y,x). >> >> There are several other possibilities for not obtaining what you >> expected. >> >> 2. Was your sequence defined over the symmetric bipolar time interval >> >> tsb = -tmax : dt : tmax % dt = 1/Fs, N is odd >> >> = dt* [ -(N-1)/2 : (N-1)/2 ]; >> >> or the asymmetric bipolar time interval >> >> tab = -(tmax+dt) : dt : tmax % N is even >> >> = dt* [ -N/2 : N/2 - 1 ]; >> >> 3. Did you use fft(xb), fft(fftshift(xb) or fft(ifftshift(xb)) ? >> All will yield the correct amplitude. However >> >> a. X = fft(ifftshift(xb)) will always yield the correct phase >> b. fft(fftshift(xb) will yield the correct phase if N is even >> and you used xb(tab) >> c. fft(xb) is not guaranteed to ever yield the correct phase. >> >> 4. Even if you used a, roundoff may have caused >> imag(X) to be very small but nonzero. >> >> An attempt at explaining how to use the correct combination >> can be obtained from my post >> >> http://groups.google.com/group/comp.soft-sys.matlab/msg/ecedcba94094f742 >> >> Hope this helps. >> >> Greg > > >> Hello Greg, > thanks, that is a good explanation. simple things get complicated. > I follow your answer. i am still unsure about why things work the way they > do.... > I feel like it is hard to trust matlab even on simple function like fft. > I am interested in looking a the phase spectrum and its changes... > Do you think I should always: > > if b is a 1D sequence, B=fft(b) and atan(imag(B)./ real(B)). his would give > the correct phase. > > For a 2D matrix a I can use A=fft2(ifftshift(a)). To plot the correct phase > using atan2(imag(A), real(A))..... > > I have tried A=fft2(a) and b=ifft2(A). it turns out that b is kind of equal > to a. > If I took the A=fft(a) and B=fft(b), they would have different phase > spectra, with strange spikes in B... > What preventive things can I do to avoid these issues? Just use > fft2(ifftshift(a)) ? > > the magnitude > > As nearly as I can tell, when you compute the DFT, the imag parts appear to be zero but with what precision? When you compute imag(A)/real(A), the results are x10^-15 and aren't zero for some reason probably having to do with the "divide" algorithm and the precision. Might that have anything to do with what you're seeing? A rotation of pi is only a flip of the sign of the sinusoid .... and +pi or -pi is the same thing. Fred
From: Greg Heath on 26 Jun 2010 00:43 On Jun 25, 11:47 am, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> wrote: > >On Jun 24, 5:12 pm, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> > >wrote: > >> Hello Forum, > > >> the fft of a rectangular, even, symmetric sequence show a nonzero real > part > >> and a zero imaginary part. All good so far. > > >> However the phase spectrum is nonzero: zero at some frequencies, pi at > some > >> frequencies , -pi at some others.... > > >> The phase spectrum is given by atan(Imag(fft)/real(fft)), > > >Incorrect. Look at the source code: Enter into the command line. > > >type angle > >type atan2 > > >> so it should give zero.... > > >> why not? > > >1. Choose your own examples and compare atan(y/x) with atan2(y,x). > > >There are several other possibilities for not obtaining what you > >expected. > > >2. Was your sequence defined over the symmetric bipolar time interval > > >tsb = -tmax : dt : tmax % dt = 1/Fs, N is odd > > > = dt* [ -(N-1)/2 : (N-1)/2 ]; > > >or the asymmetric bipolar time interval > > >tab = -(tmax+dt) : dt : tmax % N is even > > > = dt* [ -N/2 : N/2 - 1 ]; > > >3. Did you use fft(xb), fft(fftshift(xb) or fft(ifftshift(xb)) ? > > All will yield the correct amplitude. However > > >a. X = fft(ifftshift(xb)) will always yield the correct phase > >b. fft(fftshift(xb) will yield the correct phase if N is even > > and you used xb(tab) > >c. fft(xb) is not guaranteed to ever yield the correct phase. > > >4. Even if you used a, roundoff may have caused > > imag(X) to be very small but nonzero. > > >An attempt at explaining how to use the correct combination > >can be obtained from my post > > >http://groups.google.com/group/comp.soft-sys.matlab/msg/ecedcba94094f742 > > >Hope this helps. > > >Greg > >Hello Greg, > > thanks, that is a good explanation. simple things get complicated. > I follow your answer. i am still unsure about why things work the way they > do.... > I feel like it is hard to trust matlab even on simple function like fft. > I am interested in looking a the phase spectrum and its changes... > Do you think I should always: > > if b is a 1D sequence, B=fft(b) and atan(imag(B)./ real(B)). his would give > the correct phase. No. No. No. I guess you really didn't follow my answer: The use of atan is incorrect. fft phase is determined by atan2 and can be obtained from either phase = angle(B) or phase = atan2(imag(B),real(B)) Enter into the command line and compare the differences doc atan help atan doc atan2 help atan2 Now, never use atan again! > For a 2D matrix a I can use A=fft2(ifftshift(a)). provided a is defined over the "zero-centered" M X N rectangle representing yb0 = dy*[ -ceil(M-1)/2 : floor(M-1)/2 ] xb0 = dx*[ -ceil(N-1)/2 : floor(N-1)/2 ] > To plot the correct phase > using atan2(imag(A), real(A))..... % Start code close all, clear all, clc, k=0 N = 15, M = 12 dx = 0.25, dy = 0.3, xb0 = dx*[ -ceil(N-1)/2 : floor(N-1)/2 ]; yb0 = dy*[ -ceil(M-1)/2 : floor(M-1)/2 ]; xb = repmat(xb0,M,1); yb = repmat(flipud(yb0'),1,N); zb = exp(-(xb.^2 + yb.^2)); minmaxzb = minmax(zb(:)') % [ 0.0030733 0.97775] k=k+1,figure(k) contourf(xb,yb,zb) dkx = 1/(N*dx), dky = 1/(M*dy) kxb0 = dkx*[ -ceil(N-1)/2 : floor(N-1)/2 ]; kyb0 = dky*[ -ceil(M-1)/2 : floor(M-1)/2 ]; kxb = repmat(kxb0,M,1); kyb = repmat(flipud(kyb0'),1,N); Zb = fftshift(fft2(ifftshift(zb))); rZb = real(Zb); iZb = imag(Zb); aZb = abs(Zb); pZb = angle(Zb); check1 = max(maxabs(pZb-atan2(iZb,rZb)))% 0 % end code OK. It works for 2D! .... UHOH, I'm getting the feeling that you might think that the '2' in atan2 has something to do with 2D .... it doesn't. Again, use the doc and help commands above. > I have tried A=fft2(a) and b=ifft2(A). it turns out that > b is kind of equal to a. What in the world does that mean? Are you a philosophy major? Just invert the above expression for Zb: check2 = max(maxabs(zb-fftshift(ifft2(ifftshift(Zb))))) % 3.3307e-016 > If I took the A=fft(a) and B=fft(b), they would have different phase > spectra, with strange spikes in B... As defined above, a and b are 2-D. Using fft once will only transform the columns So, I'm not sure what you are getting at. > What preventive things can I do to avoid these issues? Just use > fft2(ifftshift(a)) ? No. You have to understand 1. the difference between using ifftshift and fftshift on 1-D sequences 2. the difference between using ifftshift and fftshift on 2-D sequences 3. the equivalence of fft2 and using fft twice 4. The uselessness of trying to understand phase when imag(Zb) is on the order of roundoff error. For each of the functions mentioned above, 1. Read doc ... 2. Read help ... 3. Make up a simple example and verify your understanding 4. When asking a question it sometimes helps to include the results of your simple example. Hope this helps. Greg
From: Ron N. on 26 Jun 2010 04:55 On Jun 24, 2:12 pm, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> wrote: > the fft of a rectangular, even, symmetric sequence show a nonzero real part > and a zero imaginary part. All good so far. > > However the phase spectrum is nonzero: zero at some frequencies, pi at some > frequencies , -pi at some others.... > > The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should give > zero.... Assuming the OP isn't an AI bot, it looks like I need to add another misconception to my FFT misconceptions page, that the phase of an FFT result near the size of the numerical errors in an FFT computation is relevant. IMHO. YMMV. -- rhn A.T nicholson d.0.t C-o-M http://www.nicholson.com/fftmisconceptions.html http://www.nicholson.com/rhn/dsp.html
From: Dirk Bell on 26 Jun 2010 12:43
On Jun 24, 5:12 pm, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> wrote: > Hello Forum, > > the fft of a rectangular, even, symmetric sequence show a nonzero real part > and a zero imaginary part. All good so far. > > However the phase spectrum is nonzero: zero at some frequencies, pi at some > frequencies , -pi at some others.... > > The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should give > zero.... > > why not? > > thanks > fisico32 fisico32, Your result is not surpising. If the amplitude of the spectrum is strictly real, then if the amplitude (not the magnitude) goes negative you will get a phase of +-pi (the sign depends on your software). Now consider when the ideal result should be strictly real, but there is a very small imaginary part (0 with computational error). When the real part of the of the FFT output (not the magnitude) goes negative, the sign of the computational error for the "zero" imaginary part will determine whether your software gives you a phase of essentially +pi or -pi ( when the imaginary part is <0) or essentially 0 (when the imaginary part is >0). As mentioned in a previous post, if the ideal real part is 0 and the ideal imaginary part is 0, with computational noise !=0, the result is not predictable or meaningful. Dirk |