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
nonzero phase spectrum from FFT of even function with zero imag part
Started by ●June 24, 2010
Reply by ●June 24, 20102010-06-24
On Jun 24, 5:12�pm, "fisico32" <marcoscipioni1@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? >you didn't sway the two halves with fftshift(). r b-j
Reply by ●June 24, 20102010-06-24
>On Jun 24, 5:12=A0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com> >wrote: >> Hello Forum, >> >> the fft of a rectangular, even, symmetric sequence show a nonzero realpa=>rt >> and a zero imaginary part. All good so far. >> >> However the phase spectrum is nonzero: zero at some frequencies, pi atso=>me >> frequencies , -pi at some others.... >> >> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it shouldgi=>ve >> zero.... >> >> why not? >> > >you didn't sway the two halves with fftshift(). > >r b-j >Actually, on page 101 of Understanding digital signal processing by Richard Lyon , it is said (at the bottom) that " the particular pattern of pi and -pi values is an artifact of the software used to generate that figure. A different software package may show a different pattern.... as long as the nonzero phase samples are either +pi or -pi the phase results will be correct" Is this related to the fftshift you are mentioning?
Reply by ●June 24, 20102010-06-24
On Jun 24, 2:12�pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com> wrote: ...> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should give > zero.... > > why not? > > thanks > fisico32When I try:>> x = fft([1 1 1 1 1 1])x = 6 0 0 0 0 0>> atan2(imag(x),real(x))ans = 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
Reply by ●June 24, 20102010-06-24
>On Jun 24, 2:12=A0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com> >wrote: >... >> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it shouldgi=>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...
Reply by ●June 24, 20102010-06-24
On Jun 24, 3:15�pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com> wrote:> >On Jun 24, 2:12=A0pm, "fisico32" <marcoscipioni1@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...That's why people learn to use atan2 in so many languages, including Matlab, as I did in my Matlab example. Dale B. Dalrymple
Reply by ●June 25, 20102010-06-25
On Jun 24, 5:12 pm, "fisico32" <marcoscipioni1@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
Reply by ●June 25, 20102010-06-25
On Jun 24, 6:15�pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com> wrote:> >On Jun 24, 2:12=A0pm, "fisico32" <marcoscipioni1@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
Reply by ●June 25, 20102010-06-25
>On Jun 24, 5:12 pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com> >wrote: >> Hello Forum, >> >> the fft of a rectangular, even, symmetric sequence show a nonzero realpart>> and a zero imaginary part. All good so far. >> >> However the phase spectrum is nonzero: zero at some frequencies, pi atsome>> 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
Reply by ●June 25, 20102010-06-25
fisico32 wrote:>> On Jun 24, 5:12 pm, "fisico32" <marcoscipioni1@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
