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

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 real pa=
rt
> and a zero imaginary part. All good so far.
>
> However the phase spectrum is nonzero: zero at some frequencies, pi at so=
me
> frequencies , -pi at some others....
>
> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should gi=
ve
> zero....
>
> why not?
>

you didn't sway the two halves with fftshift().

r b-j

>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 real
pa=
>rt
>> and a zero imaginary part. All good so far.
>>
>> However the phase spectrum is nonzero: zero at some frequencies, pi at
so=
>me
>> frequencies , -pi at some others....
>>
>> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should
gi=
>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?

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

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

On Jun 24, 3:15=A0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
wrote:
> >On Jun 24, 2:12=3DA0pm, "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=3D
> >ve
> >> zero....
>
> >> why not?
>
> >> thanks
> >> fisico32
>
> >When I try:
> >>> x =3D3D fft([1 1 1 1 1 1])
>
> >x =3D3D
>
> > =A0 =A0 6 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
>
> >>> atan2(imag(x),real(x))
>
> >ans =3D3D
>
> > =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 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=3D[0 1 2 3 3 2 1]; (even sequence).
> A=3Dfft(f);
> B=3Dimag(A); (it is all zeros).
> C=3Dangle(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

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

On Jun 24, 6:15=A0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
wrote:
> >On Jun 24, 2:12=3DA0pm, "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=3D
> >ve
> >> zero....
>
> >> why not?
>
> >> thanks
> >> fisico32
>
> >When I try:
> >>> x =3D3D fft([1 1 1 1 1 1])
>
> >x =3D3D
>
> > =A0 =A0 6 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0
>
> >>> atan2(imag(x),real(x))
>
> >ans =3D3D
>
> > =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 0 =A0 =A0 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=3D[0 1 2 3 3 2 1]; (even sequence).
> A=3Dfft(f);
> B=3Dimag(A); (it is all zeros).
> C=3Dangle(A) (it shows zero and pi values...)
>
> The function angle should just implement the inverse tangent...- Hide quo=
ted 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

>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 

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

On Jun 25, 11:47 am, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
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.

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

On Jun 24, 2:12=A0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
wrote:
> the fft of a rectangular, even, symmetric sequence show a nonzero real pa=
rt
> and a zero imaginary part. All good so far.
>
> However the phase spectrum is nonzero: zero at some frequencies, pi at so=
me
> frequencies , -pi at some others....
>
> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should gi=
ve
> 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

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 real pa=
rt
> and a zero imaginary part. All good so far.
>
> However the phase spectrum is nonzero: zero at some frequencies, pi at so=
me
> frequencies , -pi at some others....
>
> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should gi=
ve
> 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 !=3D0, the result is
not predictable or meaningful.

Dirk

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 real pa=
rt
> and a zero imaginary part. All good so far.
>
> However the phase spectrum is nonzero: zero at some frequencies, pi at so=
me
> frequencies , -pi at some others....
>
> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should gi=
ve
> 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 !=3D0, the result
is
not predictable or meaningful.


Dirk

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 real pa=
rt
> and a zero imaginary part. All good so far.
>
> However the phase spectrum is nonzero: zero at some frequencies, pi at so=
me
> frequencies , -pi at some others....
>
> The phase spectrum is given by atan(Imag(fft)/real(fft)), so it should gi=
ve
> 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).

Dirk

On Jun 24, 5:48=A0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
wrote:
> >On Jun 24, 5:12=3DA0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com=
>
....
> >> why not?
>

r b-j said (with spelling correction)

> >you didn't swap the two halves with fftshift().
>
> >r b-j
>
> Actually, on page 101 of Understanding digital signal processing by Richa=
rd
> 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 th=
e
> nonzero phase samples are either +pi or -pi the phase results will be
> correct"
>
> Is this related to the fftshift you are mentioning?

i dunno if Rick's thing is related until i dig the book out of the
correct box (since moving it from Waltham Mass back to Vermont, all of
my books are in boxes).

if Rick is saying what i think you are crediting him for saying, then
i think that Rick is sorta incorrect.  this *is* directly related to
failing to use fftshift() (or it's counterpart in whatever language or
package you are using) before applying the FFT.  if you have a vector
of data (with indices 0 to N-1) that is extracted from some continuous
stream of data much longer than N, and if you apply a window (like the
Hann or Hamming or Kaiser or whatever) to that vector such that the
middle of that vector (around index N/2) is where the middle and
maximum of the window is, and send that directly into the FFT without
swapping the first and last halves, then if you examine the phase, you
will see a phase shift of +/- pi (doesn't matter which, since they are
2*pi apart) added to *only* the odd numbered values of X[k].

it's not just a display thing.  it is also an analysis thing.  if
phase continuity is important (and if you're bothering to unwrap
phase, continuity *is* important) then it is necessary to swap the
first and last halves of the data before it goes into the FFT.  that's
how you avoid that icky problem.

r b-j

On Jun 26, 9:54 pm, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Jun 24, 5:48 pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
> wrote:
>
> > >On Jun 24, 5:12=A0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
> ...
> > >> why not?
>
> r b-j said (with spelling correction)
>
> > >you didn't swap 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?

No, it is not.

In addition, it is not related to your problem either.

> i dunno if Rick's thing is related until i dig the book out of the
> correct box (since moving it from Waltham Mass back to Vermont, all of
> my books are in boxes).
>
> if Rick is saying what i think you are crediting him for saying, then
> i think that Rick is sorta incorrect.  this *is* directly related to
> failing to use fftshift() (or it's counterpart in whatever language or
> package you are using) before applying the FFT.

No, it is not.

> if you have a vector
> of data (with indices 0 to N-1) that is extracted from some continuous
> stream of data much longer than N, and if you apply a window (like the
> Hann or Hamming or Kaiser or whatever) to that vector such that the
> middle of that vector (around index N/2) is where the middle and
> maximum of the window is, and send that directly into the FFT without
> swapping the first and last halves, then if you examine the phase, you
> will see a phase shift of +/- pi (doesn't matter which, since they are
> 2*pi apart) added to *only* the odd numbered values of X[k].
>
> it's not just a display thing.  it is also an analysis thing.  if
> phase continuity is important (and if you're bothering to unwrap
> phase, continuity *is* important) then it is necessary to swap the
> first and last halves of the data before it goes into the FFT.  that's
> how you avoid that icky problem.

For the OP's particular problem, all of this is irrelevant.

His posted example was

clear all,  clc

f = [0 1 2 3 3 2 1]' ; F = fft(f); realF = real(F); imagF = imag(F);
phaseF = angle(F); wrongphase = atan(imagF./realF);
rightphase = atan2(imagF,realF);

summary1 = [f F realF imagF phaseF wrongphase rightphase]

% f      F              realF     imagF  phaseF  wrongphase rightphase
% 0     12             12          0           0
0                0
% 1   -5.0489    -5.0489      0       3.1416          0
3.1416
% 2   -0.30798  -0.30798    0       3.1416          0           3.1416
% 3   -0.6431    -0.6431      0       3.1416         0
3.1416
% 3   -0.6431    -0.6431      0       3.1416         0
3.1416
% 2   -0.30798  -0.30798     0      3.1416          0           3.1416
% 1   -5.0489    -5.0489     0        3.1416         0
3.1416

disp('Even if the imaginary component is just floating point error:')

F = F + i*eps*rand(7,1); realF = real(F); imagF = imag(F);
phaseF = angle(F); wrongphase = atan(imagF./realF);
rightphase = atan2(imagF,realF);

summary2 = [realF imagF phaseF wrongphase rightphase]

% realF         imagF       phaseF        wrongphase    rightphase
% 12         2.1097e-016   1.7581e-017   1.7581e-017
1.7581e-017
% -5.0489    5.1323e-017   3.1416       -1.0165e-017    3.1416
% -0.30798   1.3475e-016   3.1416       -4.3752e-016
3.1416
% -0.6431    1.0791e-016   3.1416       -1.678e-016     3.1416
% -0.6431    1.9791e-016   3.1416       -3.0774e-016    3.1416
% -0.30798   1.6922e-016   3.1416       -5.4945e-016    3.1416
% -5.0489    1.0136e-016   3.1416       -2.0075e-017
3.1416

disp(['Therefore the OPs problem is not recognizing that phase is';
      ' the result of the atan2 fucntion, not the atan function !'])

Hope this helps.

Greg

On Jun 26, 4:55=A0am, "Ron N." <rhnlo...@yahoo.com> wrote:
> On Jun 24, 2:12=A0pm, "fisico32" <marcoscipioni1@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.

Also add the misconception that caused this thread:

fft phase is obtained using atan instead of atan2.

Hope this helps.

Greg

On Jun 27, 11:41=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 26, 9:54 pm, robert bristow-johnson <r...@audioimagination.com>
> wrote:
>
>
>
> > On Jun 24, 5:48 pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
> > wrote:
>
> > > >On Jun 24, 5:12=3DA0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail=
..com>
> > ...
> > > >> why not?
>
> > r b-j said (with spelling correction)
>
> > > >you didn't swap the two halves with fftshift().
>
> > > >r b-j
>
> > > Actually, on page 101 of Understanding digital signal processing by R=
ichard
> > > 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 figur=
e. A
> > > different software package may show a different pattern.... as long a=
s the
> > > nonzero phase samples are either +pi or -pi the phase results will be
> > > correct"
>
> > > Is this related to the fftshift you are mentioning?
>
> No, it is not.
>
> In addition, it is not related to your problem either.
>
> > i dunno if Rick's thing is related until i dig the book out of the
> > correct box (since moving it from Waltham Mass back to Vermont, all of
> > my books are in boxes).
>
> > if Rick is saying what i think you are crediting him for saying, then
> > i think that Rick is sorta incorrect. =A0this *is* directly related to
> > failing to use fftshift() (or it's counterpart in whatever language or
> > package you are using) before applying the FFT.
>
> No, it is not.
>

so Greg, just to avoid goofy things that happens when the sequence is
odd-lengthed, let's make the sequence a little longer and make it even-
lengthed (and keep the symmetry).


so can you tell me what you might expect with this sequence?:


    x =3D [0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1]';

    X =3D fft(x);

    plot(angle(X + 1e-10));


versus this one:


    x =3D [8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7]';

    X =3D fft(x);

    plot(angle(X + 1e-10));


what might you expect to see happen to the phases of the odd-indexed
elements of X (actually, since it's stupid MATLAB that puts DC into
X(1), the question is what happens to the phases of the even-indexed
elements of X)?

the OP may have problems with the angle of 1+j0 vs. -1+j0 and MATLAB
(or any numerical package) might have a problem with angle(-0 + j0)
vs. angle(+0 + j0) which is why i added 1e-10 to the real part of X.


> Hope this helps.

it doesn't.  but there's always hope.

also, consider posting column-aligned tables with a mono-spaced font.
you cannot assume that the proportional width font you're using is the
same that anyone else does.  i didn't even bother to look at it, i
wasn't going to try to align the columns.

r b-j

On Jun 27, 2:00 pm, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Jun 27, 11:41 am, Greg Heath <he...@alumni.brown.edu> wrote:
> > On Jun 26, 9:54 pm, robert bristow-johnson <r...@audioimagination.com>
> > wrote:
> > > On Jun 24, 5:48 pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
> > > wrote:
> > > > >On Jun 24, 5:12=A0pm, "fisico32" <marcoscipioni1@n_o_s_p_a_m.gmail.com>
> > > ...
> > > > >> why not?
> > > r b-j said (with spelling correction)
> > > > >you didn't swap the two halves with fftshift().

Original OP (fisico32) question: if F = fft(f) and
imag(F) = 0, why isn't angle(F) = 0 as predicted by
atan?

Answer: angle is calculated using atan2 in the half
closed range (-pi,+pi], not by using atan which only
has the closed range [-pi/2,+pi/2].

My example based on the original question:

eps       = 2.2204e-016  % (MATLAB machine epsilon)
correct   = atan2( [0 0 0 0 0]',[-1 -eps 0 eps 1]');
incorrect = atan( [0 0 0 0 0]'./[-1 -eps 0 eps 1]');

"correct"         "incorrect"(sometimes)
3.1416                     0   % incorrect
3.1416                     0   % incorrect
0                          NaN  % incorrect (0/0)
0                            0
0                            0

Therefore, imag(F) = 0, real(F) <0 ==> angle(X) = pi.

A little OT excursion for a moment: Assume real(F) = 0

correct   = atan2([-1 -eps 0 eps 1]' ,  [0 0 0 0 0]')
incorrect = atan( [-1 -eps 0 eps 1]' ./ [0 0 0 0 0]')

"correct"          "incorrect"(sometimes)
 -1.5708                   -1.5708
 -1.5708                   -1.5708
       0                          NaN    % incorrect (0/0)
  1.5708                    1.5708
  1.5708                    1.5708

Now, back to imag(F) = 0:

Since the OP did not give an example and there are
other causes of unexpected nonzero phase, there
were a variety of possibilities mentioned in replies:

1. misuse of fftshift
2. atan
3. floating point roundoff error
4. original array misalignment(not caused by shifting)

Later, the OP gave an example that clearly indicated
his confusion was caused by the atan (instead of atan2)
misconception.

This does not mean that I think the other causes
don't occur. They do. I just wanted to emphasize that
the other causes are not the answer to the example
posted by the OP.

Now, let's talk fftshift:

> so Greg, just to avoid goofy things that happens when the sequence is
> odd-lengthed, let's make the sequence a little longer and make it even-
> lengthed (and keep the symmetry).
>
> so can you tell me what you might expect with this sequence?:

Indices 1 and 2 added for clarity (GEH)
The 1e-10 fudge is unnecessary with MATLAB.

>     x1 = [0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1]';
>
>     X1 = fft(x1);
>     plot(angle(X1 + 1e-10));
>
> versus this one:
>
>     x2 = [8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7]';
>
>     X2 = fft(x2);

>     plot(angle(X2 + 1e-10));

Since t is not explicitly defined it must be
assumed that

dt = 1/Fs, t = dt*[0:N-1];

Both can be considered full N = 16 point periods
of the same periodic waveform except for a mutual
half period delay of t0 = (N/2)*dt seconds.

Both can be considered sampled from even periodic
functions so imag(X2) = imag(X1) = 0.

Furthermore, if

f = (Fs/N)*[0:N-1];

angle(X2)-angle(X1) = 2*pi*f*t0 = pi*[0:N-1}

Therefore, the real parts will differ in
sign every other point.

> what might you expect to see happen to the phases of the odd-indexed
> elements of X (actually, since it's stupid MATLAB that puts DC into
> X(1), the question is what happens to the phases of the even-indexed
> elements of X)?

real(X) >= 0 ==> angle(X) = 0
real(X) < 0  ==> angle(X) = pi

However,  what property of x guarantees that

real(X) >= 0 ?

> the OP may have problems with the angle of 1+j0 vs. -1+j0

Only if atan is used.

>and MATLAB
> (or any numerical package) might have a problem with angle(-0 + j0)
> vs. angle(+0 + j0) which is why i added 1e-10 to the real part of X.

No problems with MATLAB.

> > Hope this helps.
>
> it doesn't.  but there's always hope.

Hopefully, my preface in this reply made it clear that
the atan misconception was the cause of nonzero phase
in the example posted by the OP. However, as listed
there are other possible causes for other examples.

> also, consider posting column-aligned tables with a mono-spaced font.
> you cannot assume that the proportional width font you're using is the
> same that anyone else does.  i didn't even bother to look at it, i
> wasn't going to try to align the columns.

I post in Google groups. AFAIK there is no way do do that.
However, I will check.

Hope this helps.

Greg

On Jun 28, 1:12=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 27, 2:00 pm, robert bristow-johnson <r...@audioimagination.com>
> wrote:
>

>
> Original OP (fisico32) question: if F =3D fft(f) and
> imag(F) =3D 0, why isn't angle(F) =3D 0 as predicted by
> atan?
>
> Answer: angle is calculated using atan2 in the half
> closed range (-pi,+pi], not by using atan which only
> has the closed range [-pi/2,+pi/2].
>

> Since the OP did not give an example and there are
> other causes of unexpected nonzero phase, there
> were a variety of possibilities mentioned in replies:
>
> 1. misuse of fftshift

what would that be?  (though i dunno how i would use it for N odd.)

> 2. atan

we both agree with that.  and atan() would not ever put an angle
outside of quadrants 1 and 4.

> 3. floating point roundoff error

that's the +/- 0 + j0 issue, which is why i added 1e-10.

> 4. original array misalignment(not caused by shifting)

well that misalignment is what i was addressing.  and how is alignment
mis-done without shifting?


> Now, let's talk fftshift:
>
> > so Greg, just to avoid goofy things that happens when the sequence is
> > odd-lengthed, let's make the sequence a little longer and make it even-
> > lengthed (and keep the symmetry).
>
> > so can you tell me what you might expect with this sequence?:
>
> Indices 1 and 2 added for clarity (GEH)

good idea.

> The 1e-10 fudge is unnecessary with MATLAB.

it was needed for Octave (which i run on my laptop).  haven't checked
it with MATLAB yet.

> > =A0 =A0 x1 =3D [0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1]';
> > =A0 =A0 X1 =3D fft(x1);
> > =A0 =A0 plot(angle(X1 + 1e-10));
>
> > versus this one:
>
> > =A0 =A0 x2 =3D [8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7]';
> > =A0 =A0 X2 =3D fft(x2);
> > =A0 =A0 plot(angle(X2 + 1e-10));


angle(X2 + 1e-10) should have been a straight horizontal line at 0.

angle(X1 + 1e-10) would be equal to +/- pi for only even indices (in
MATLAB/Octave) or for only odd indices for normal people (who do not
excuse MATLAB for screwing their indices up).

one has nice continuous phase (because it is centered about t=3D0) and
the other is delayed by N/2 samples with the associated phase shift
involved.  since the angle() continues to be wrapped, it doesn't
advance by -pi for each sample, but wraps instead.

> Since t is not explicitly defined it must be
> assumed that
>
> dt =3D 1/Fs, t =3D dt*[0:N-1];

sure, but i wasn't worrying about continuous time.  i was just leaving
it in discrete time.  however, i *do* worry about the inherent
periodicity of the DFT.

> Both can be considered full N =3D 16 point periods
> of the same periodic waveform except for a mutual
> half period delay of t0 =3D (N/2)*dt seconds.
>
> Both can be considered sampled from even periodic
> functions so imag(X2) =3D imag(X1) =3D 0.
>
> Furthermore, if
>
> f =3D (Fs/N)*[0:N-1];
>
> angle(X2)-angle(X1) =3D 2*pi*f*t0 =3D pi*[0:N-1}
>
> Therefore, the real parts will differ in
> sign every other point.

yah.  and with the imaginary part always zero (as was the case for the
OP) what does that do to the phase when you do the proper atan2()
thingie?

> > what might you expect to see happen to the phases of the odd-indexed
> > elements of X (actually, since it's stupid MATLAB that puts DC into
> > X(1), the question is what happens to the phases of the even-indexed
> > elements of X)?
>
> real(X) >=3D 0 =3D=3D> angle(X) =3D 0

for the even indices (or odd indices in MATLABland).

> real(X) < 0 =A0=3D=3D> angle(X) =3D pi

for the odd indices (or even inside of MATLAB).

> However, =A0what property of x guarantees that
>
> real(X) >=3D 0 ?

none in the example that i know of.  it oscillates every other sample.

> > the OP may have problems with the angle of 1+j0 vs. -1+j0
>
> Only if atan is used.

yup, we agree about that.

> > and MATLAB
> > (or any numerical package) might have a problem with angle(-0 + j0)
> > vs. angle(+0 + j0) which is why i added 1e-10 to the real part of X.
>
> No problems with MATLAB.

well, it was with my Octave implementation.  (it wasn't really +/- 0,
which IEEE 754 defines, but more like +/- 10^-15 from floating-point
errors.  it was enough to flip the angle back and forth.)


> Hopefully, my preface in this reply made it clear that
> the atan misconception was the cause of nonzero phase
> in the example posted by the OP.

but it's not the only problem.  and the other "misalignment" issue is
about the only reason i ever use fftshift().

> However, as listed
> there are other possible causes for other examples.
>
> > also, consider posting column-aligned tables with a mono-spaced font.
> > you cannot assume that the proportional width font you're using is the
> > same that anyone else does. =A0i didn't even bother to look at it, i
> > wasn't going to try to align the columns.
>
> I post in Google groups. AFAIK there is no way do do that.
> However, I will check.

click the options for the thread.  then hit "Fixed font".

L8r,

r b-j

On Jun 28, 2:41=A0pm, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Jun 28, 1:12=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> > On Jun 27, 2:00 pm, robert bristow-johnson <r...@audioimagination.com>
> > wrote:
>
> > Original OP (fisico32) question: if F =3D fft(f) and
> > imag(F) =3D 0, why isn't angle(F) =3D 0 as predicted by
> > atan?
>
> > Answer: angle is calculated using atan2 in the half
> > closed range (-pi,+pi], not by using atan which only
> > has the closed range [-pi/2,+pi/2].
>
> > Since the OP did not give an example and there are
> > other causes of unexpected nonzero phase, there
> > were a variety of possibilities mentioned in replies:
>
> > 1. misuse of fftshift
>
> what would that be? =A0(though i dunno how i would use it for N odd.)

For N odd, using fftshift instead of fftshift and vice versa.

The subscript b refers to zero centered "b"ipolar
time and frequency intervals proportional
(via dt=3D1/Fs or df=3DFs/N) to

[-ceil((N-1)/2) : floor((N-1)/2)]

x =3D ifftshift(xb)
X =3D fft(x)
Xb =3D fftshift(X)

X =3D ifftshift(Xb)
x =3D ifft(X)
xb =3D fftshift(x)


> > 2. atan
>
> we both agree with that. =A0and atan() would not ever put an angle
> outside of quadrants 1 and 4.
>
> > 3. floating point roundoff error
>
> that's the +/- 0 + j0 issue, which is why i added 1e-10.

The problem in MATLAB is not when real(F) or imag(F)
is exactly zero, but when it is all bipolar roundoff error.

The only solution that I have found for that is to
threshold based on abs(imag(F))./max(abs(F))

> > 4. original array misalignment(not caused by shifting)
>
> well that misalignment is what i was addressing. =A0and how is alignment
> mis-done without shifting?

I meant not caused by fftshift or ifftshift but
by coding error or ignorance. For example,
defining

t =3D [0:dt:T]  instead of t =3D[0:dt:T-dt]
or
tb =3D [-T/2:dt:T/2] instead of tb =3D [-T/2:dt:T/2-dt]

> > Now, let's talk fftshift:
>
> > > so Greg, just to avoid goofy things that happens when the sequence is
> > > odd-lengthed, let's make the sequence a little longer and make it eve=
n-
> > > lengthed (and keep the symmetry).
>
> > > so can you tell me what you might expect with this sequence?:
>
> > Indices 1 and 2 added for clarity (GEH)
>
> good idea.
>
> > The 1e-10 fudge is unnecessary with MATLAB.
>
> it was needed for Octave (which i run on my laptop). =A0haven't checked
> it with MATLAB yet.
>
> > > =A0 =A0 x1 =3D [0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1]';
> > > =A0 =A0 X1 =3D fft(x1);
> > > =A0 =A0 plot(angle(X1 + 1e-10));
>
> > > versus this one:
>
> > > =A0 =A0 x2 =3D [8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7]';
> > > =A0 =A0 X2 =3D fft(x2);
> > > =A0 =A0 plot(angle(X2 + 1e-10));
>
> angle(X2 + 1e-10) should have been a straight horizontal line at 0.
OK
> angle(X1 + 1e-10) would be equal to +/- pi for only even indices (in
> MATLAB/Octave) or for only odd indices for normal people (who do not
> excuse MATLAB for screwing their indices up).

I got 0 for odd indices and +pi for odd.

However, if I had gotten bipolar roundoff instead
of exactly 0 for imag(X), I would probably have
gotten some -pi.

> one has nice continuous phase (because it is centered about t=3D0) and
> the other is delayed by N/2 samples with the associated phase shift
> involved. =A0since the angle() continues to be wrapped, it doesn't
> advance by -pi for each sample, but wraps instead.

I didn't unwrap the angle.

> > Since t is not explicitly defined it must be
> > assumed that
>
> > dt =3D 1/Fs, t =3D dt*[0:N-1];
>
> sure, but i wasn't worrying about continuous time. =A0i was just
> leaving it in discrete time. =A0

That is discrete time: t =3D [0,dt,2*dt,...(N-1)*dt];

however, i *do* worry about the inherent
> periodicity of the DFT.
>
> > Both can be considered full N =3D 16 point periods
> > of the same periodic waveform except for a mutual
> > half period delay of t0 =3D (N/2)*dt seconds.
>
> > Both can be considered sampled from even periodic
> > functions so imag(X2) =3D imag(X1) =3D 0.
>
> > Furthermore, if
>
> > f =3D (Fs/N)*[0:N-1];
>
> > angle(X2)-angle(X1) =3D 2*pi*f*t0 =3D pi*[0:N-1}
>
> > Therefore, the real parts will differ in
> > sign every other point.
>
> yah. =A0and with the imaginary part always zero (as was the case for the
> OP) what does that do to the phase when you do the proper atan2()
> thingie?
>
> > > what might you expect to see happen to the phases of the odd-indexed
> > > elements of X (actually, since it's stupid MATLAB that puts DC into
> > > X(1), the question is what happens to the phases of the even-indexed
> > > elements of X)?
>
> > real(X) >=3D 0 =3D=3D> angle(X) =3D 0
>
> for the even indices (or odd indices in MATLABland).
>
> > real(X) < 0 =A0=3D=3D> angle(X) =3D pi
>
> for the odd indices (or even inside of MATLAB).
>
> > However, =A0what property of x guarantees that
>
> > real(X) >=3D 0 ?
>
> none in the example that i know of. =A0it oscillates every other sample.
>
> > > the OP may have problems with the angle of 1+j0 vs. -1+j0
>
> > Only if atan is used.
>
> yup, we agree about that.
>
> > > and MATLAB
> > > (or any numerical package) might have a problem with angle(-0 + j0)
> > > vs. angle(+0 + j0) which is why i added 1e-10 to the real part of X.
>
> > No problems with MATLAB.
>
> well, it was with my Octave implementation. =A0(it wasn't really +/- 0,
> which IEEE 754 defines, but more like +/- 10^-15 from floating-point
> errors. =A0it was enough to flip the angle back and forth.)
>
> > Hopefully, my preface in this reply made it clear that
> > the atan misconception was the cause of nonzero phase
> > in the example posted by the OP.
>
> but it's not the only problem. =A0and the other "misalignment" issue is
> about the only reason i ever use fftshift().
>
> > However, as listed
> > there are other possible causes for other examples.
>
> > > also, consider posting column-aligned tables with a mono-spaced font.
> > > you cannot assume that the proportional width font you're using is th=
e
> > > same that anyone else does. =A0i didn't even bother to look at it, i
> > > wasn't going to try to align the columns.
>
> > I post in Google groups. AFAIK there is no way do do that.
> > However, I will check.
>
> click the options for the thread. =A0then hit "Fixed font".

Duh. Well, you can't expect me to have known that.
I've only been using Google Groups since 1996.

Right?

Greg

On Jun 28, 4:46=A0pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 28, 2:41=A0pm, robert bristow-johnson <r...@audioimagination.com>
> wrote:
-----SNIP
> > > > also, consider posting column-aligned tables with a mono-spaced fon=
t.
> > > > you cannot assume that the proportional width font you're using is =
the
> > > > same that anyone else does. =A0i didn't even bother to look at it, =
i
> > > > wasn't going to try to align the columns.
>
> > > I post in Google groups. AFAIK there is no way do do that.
> > > However, I will check.
>
> > click the options for the thread. =A0then hit "Fixed font".
>
> Duh. Well, you can't expect me to have known that.
> I've only been using Google Groups since 1996.
>
> Right?

I have posted in the past with precisely aligned columns.
It has to be one of the following:

1. Typed them into the reply window
2. Cut and pasted from a MATLAB *.m file
3. Cut and pasted from a WINDOZE *.txt file
4. Cut and pasted from a WINDOZE *.rtf file

I'll figure it out.

Greg

On Jun 28, 4:46=A0pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jun 28, 2:41=A0pm, robert bristow-johnson <r...@audioimagination.com>
> wrote:
>
>
>
> > On Jun 28, 1:12=A0am, Greg Heath <he...@alumni.brown.edu> wrote:
> > > On Jun 27, 2:00 pm, robert bristow-johnson <r...@audioimagination.com=
>
> > > wrote:
>
> > > Original OP (fisico32) question: if F =3D fft(f) and
> > > imag(F) =3D 0, why isn't angle(F) =3D 0 as predicted by
> > > atan?
>
> > > Answer: angle is calculated using atan2 in the half
> > > closed range (-pi,+pi], not by using atan which only
> > > has the closed range [-pi/2,+pi/2].
>
> > > Since the OP did not give an example and there are
> > > other causes of unexpected nonzero phase, there
> > > were a variety of possibilities mentioned in replies:
>
> > > 1. misuse of fftshift
>
> > what would that be? =A0(though i dunno how i would use it for N odd.)
>
> For N odd, using fftshift instead of fftshift and vice versa.

my problem with MATLAB is that i always use fftshift() instead of
fftshift().  but i never had N odd for the fft().

> The subscript b refers to zero centered "b"ipolar
> time and frequency intervals proportional
> (via dt=3D1/Fs or df=3DFs/N) to
>
> [-ceil((N-1)/2) : floor((N-1)/2)]
>
> x =3D ifftshift(xb)
> X =3D fft(x)
> Xb =3D fftshift(X)
>
> X =3D ifftshift(Xb)
> x =3D ifft(X)
> xb =3D fftshift(x)
>
> > > 2. atan
>
> > we both agree with that. =A0and atan() would not ever put an angle
> > outside of quadrants 1 and 4.
>
> > > 3. floating point roundoff error
>
> > that's the +/- 0 + j0 issue, which is why i added 1e-10.
>
> The problem in MATLAB is not when real(F) or imag(F)
> is exactly zero, but when it is all bipolar roundoff error.
>
> The only solution that I have found for that is to
> threshold based on abs(imag(F))./max(abs(F))
>
> > > 4. original array misalignment(not caused by shifting)
>
> > well that misalignment is what i was addressing. =A0and how is alignmen=
t
> > mis-done without shifting?
>
> I meant not caused by fftshift or ifftshift but
> by coding error or ignorance. For example,
> defining
>
> t =3D [0:dt:T] =A0instead of t =3D[0:dt:T-dt]

the N-1 issue.  hard to catch in MATLAB that likes to count from 1 to
N instead of 0 to N-1.

> or
> tb =3D [-T/2:dt:T/2] instead of tb =3D [-T/2:dt:T/2-dt]

it would be even nicer if MATLAB allowed negative indices in addition
to 0.  but in MATLAB-land, indices are only positive integers.

>
>
>
> > > Now, let's talk fftshift:
>
> > > > so Greg, just to avoid goofy things that happens when the sequence =
is
> > > > odd-lengthed, let's make the sequence a little longer and make it e=
ven-
> > > > lengthed (and keep the symmetry).
>
> > > > so can you tell me what you might expect with this sequence?:
>
> > > Indices 1 and 2 added for clarity (GEH)
>
> > good idea.
>
> > > The 1e-10 fudge is unnecessary with MATLAB.
>
> > it was needed for Octave (which i run on my laptop). =A0haven't checked
> > it with MATLAB yet.
>
> > > > =A0 =A0 x1 =3D [0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1]';
> > > > =A0 =A0 X1 =3D fft(x1);
> > > > =A0 =A0 plot(angle(X1 + 1e-10));
>
> > > > versus this one:
>
> > > > =A0 =A0 x2 =3D [8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7]';
> > > > =A0 =A0 X2 =3D fft(x2);
> > > > =A0 =A0 plot(angle(X2 + 1e-10));
>
> > angle(X2 + 1e-10) should have been a straight horizontal line at 0.
> OK
> > angle(X1 + 1e-10) would be equal to +/- pi for only even indices (in
> > MATLAB/Octave) or for only odd indices for normal people (who do not
> > excuse MATLAB for screwing their indices up).
>
> I got 0 for odd indices and +pi for odd.

if you meant "even" for the first "odd", then you ain't counting the
MATLAB way (you're counting the Dykstra way:
 http://userweb.cs.utexas.edu/users/EWD/transcriptions/EWD08xx/EWD831.html
).

> However, if I had gotten bipolar roundoff instead
> of exactly 0 for imag(X), I would probably have
> gotten some -pi.

the issue is the toggling of +/- pi vs. 0 for every other sample (that
can be avoided if you use fftshift to center things properly).  and it
isn't an issue of roundoff error around 0.

> > one has nice continuous phase (because it is centered about t=3D0) and
> > the other is delayed by N/2 samples with the associated phase shift
> > involved. =A0since the angle() continues to be wrapped, it doesn't
> > advance by -pi for each sample, but wraps instead.
>
> I didn't unwrap the angle.

neither did i, but the effect of *not* unwrapping the angle is that
every other sample (or bin) of X1 gets pi or -pi added to it.  if
either of us *did* unwrap it, you would see a linear ramp of phase
with slope of -pi per bin that would be the appropriate phase shift
for a delay of x1[] (relative to x2[]) of N/2 samples.

the whole idea of swapping the earlier and latter halves of x1[] is to
undo that delay of N/2 samples and to center the windowed piece of
data about x1[0] (which is what x2[] is).

> > > Since t is not explicitly defined it must be
> > > assumed that
>
> > > dt =3D 1/Fs, t =3D dt*[0:N-1];
>
> > sure, but i wasn't worrying about continuous time. =A0i was just
> > leaving it in discrete time. =A0
>
> That is discrete time: t =3D [0,dt,2*dt,...(N-1)*dt];

again, i really don't give a rat's ass about dt.  it's just a positive
and real scaling factor (or a measure of time or something else in
terms of an anthropocentric definition of the unit) and it doesn't
change the phase angle.  dt could just as well be 1.

> > however, i *do* worry about the inherent periodicity of the DFT.

....

>
> > > > also, consider posting column-aligned tables with a mono-spaced fon=
t.
> > > > you cannot assume that the proportional width font you're using is =
the
> > > > same that anyone else does. =A0i didn't even bother to look at it, =
i
> > > > wasn't going to try to align the columns.
>
> > > I post in Google groups. AFAIK there is no way do do that.
> > > However, I will check.
>
> > click the options for the thread. =A0then hit "Fixed font".

actually it says "Fixed text" until you hit it, then it says "Fixed
font" (but "Proportional text" remains unchanged, other than the
underline).

> Duh. Well, you can't expect me to have known that.
> I've only been using Google Groups since 1996.
>
> Right?

right.

r b-j

On 6/28/2010 11:32 PM, robert bristow-johnson wrote:
> On Jun 28, 4:46 pm, Greg Heath<he...@alumni.brown.edu>  wrote:

   ...

>> For N odd, using fftshift instead of fftshift and vice versa.
>
> my problem with MATLAB is that i always use fftshift() instead of
> fftshift().  ...

Are you guys kidding, or what?

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

On Jun 29, 3:25=A0am, Jerry Avins <j...@ieee.org> wrote:
> On 6/28/2010 11:32 PM, robert bristow-johnson wrote:
>
> > On Jun 28, 4:46 pm, Greg Heath<he...@alumni.brown.edu> =A0wrote:
>
> =A0 =A0...
>
> >> For N odd, using fftshift instead of fftshift and vice versa.
>
> > my problem with MATLAB is that i always use fftshift() instead of
> > fftshift(). =A0...
>
> Are you guys kidding, or what?

No, not both. Just one. The other is typo prone.

Greg

On Jun 28, 11:32 pm, robert bristow-johnson
<r...@audioimagination.com> wrote:
> On Jun 28, 4:46 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> > On Jun 28, 2:41 pm, robert bristow-johnson <r...@audioimagination.com>
> > wrote:
> > > On Jun 28, 1:12 am, Greg Heath <he...@alumni.brown.edu> wrote:
> > > > On Jun 27, 2:00 pm, robert bristow-johnson <r...@audioimagination.com>
> > > > wrote:
>
> > > > Original OP (fisico32) question: if F = fft(f) and
> > > > imag(F) = 0, why isn't angle(F) = 0 as predicted by
> > > > atan?
>
> > > > Answer: angle is calculated using atan2 in the half
> > > > closed range (-pi,+pi], not by using atan which only
> > > > has the closed range [-pi/2,+pi/2].
>
> > > > Since the OP did not give an example and there are
> > > > other causes of unexpected nonzero phase, there
> > > > were a variety of possibilities mentioned in replies:
>
> > > > 1. misuse of fftshift

For N odd, using fftshift instead of ifftshift and vice versa.

> > The subscript b refers to zero centered "b"ipolar
> > time and frequency intervals proportional
> > (via dt=1/Fs or df=Fs/N) to
>
> > [-ceil((N-1)/2) : floor((N-1)/2)]
>
> > x = ifftshift(xb)
> > X = fft(x)
> > Xb = fftshift(X)
>
> > X = ifftshift(Xb)
> > x = ifft(X)
> > xb = fftshift(x)
>
> > > > 2. atan
>
> > > > 3. floating point roundoff error
>
> > The problem in MATLAB is not when real(F) or imag(F)
> > is exactly zero, but when it is all bipolar roundoff error.
>
> > The only solution that I have found for that is to
> > threshold based on abs(imag(F))./max(abs(F))
>
> > > > 4. original array misalignment(not caused by shifting)
>
> > I meant not caused by fftshift or ifftshift but
> > by coding error or ignorance. For example,
> > defining
>
> > t = [0:dt:T]  instead of t =[0:dt:T-dt]
>
> the N-1 issue.  hard to catch in MATLAB that likes to count from 1 to
> N instead of 0 to N-1.
>
> > or
> > tb = [-T/2:dt:T/2] instead of tb = [-T/2:dt:T/2-dt]
>
> it would be even nicer if MATLAB allowed negative indices in addition
> to 0.  but in MATLAB-land, indices are only positive integers.
>
> > > > Now, let's talk fftshift:
>
> > > > > so Greg, just to avoid goofy things that happens when the sequence is
> > > > > odd-lengthed, let's make the sequence a little longer and make it even-
> > > > > lengthed (and keep the symmetry).
>
> > > > > so can you tell me what you might expect with this sequence?:
>
> > > > >     x1 = [0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1]';
>
> > > > > versus this one:
>
> > > > >     x2 = [8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7]';
>
> > > angle(X2 + 1e-10) should have been a straight horizontal line at 0.
> > OK
> > > angle(X1 + 1e-10) would be equal to +/- pi for only even indices (in
> > > MATLAB/Octave) or for only odd indices for normal people (who do not
> > > excuse MATLAB for screwing their indices up).

I got 0 for the odd MATLAB indices and +pi for the even.

> > However, if I had gotten bipolar roundoff instead
> > of exactly 0 for imag(X), I would probably have
> > gotten some -pi.

Because of the half closed definition -pi < atan2 <= pi,
I would probably have gotten some within roundoff error
of -pi

> the issue is the toggling of +/- pi vs. 0 for every other sample (that
> can be avoided if you use fftshift to center things properly).  and it
> isn't an issue of roundoff error around 0.

As a wiseass man named Bill Clinton once said:

....It depends on what you mean by "issue" ...

Any waveform used in fft will yield a phase that is referenced
to the coordinate location corresponding to the first element.
If the imaginary part of the transform is exactly zero, the phase
will either be 0  (for real(F) >=0) or +pi (for real(F) < 0).

If a waveform is positive symmetric, the transform can be made
to be nonnegative real by changing the phase reference to the
symmetry point.

Not sure why this warrants an "issue" label.

> > > one has nice continuous phase (because it is centered about t=0) and
> > > the other is delayed by N/2 samples with the associated phase shift
> > > involved.  since the angle() continues to be wrapped, it doesn't
> > > advance by -pi for each sample, but wraps instead.
>
> > I didn't unwrap the angle.
>
> neither did i, but the effect of *not* unwrapping the angle is that
> every other sample (or bin) of X1 gets pi or -pi added to it.  if
> either of us *did* unwrap it, you would see a linear ramp of phase
> with slope of -pi per bin that would be the appropriate phase shift
> for a delay of x1[] (relative to x2[]) of N/2 samples.
>
> the whole idea of swapping the earlier and latter halves of x1[] is to
> undo that delay of N/2 samples and to center the windowed piece of
> data about x1[0] (which is what x2[] is).
>
> > > > Since t is not explicitly defined it must be
> > > > assumed that
>
> > > > dt = 1/Fs, t = dt*[0:N-1];
>
> > > sure, but i wasn't worrying about continuous time.  i was just
> > > leaving it in discrete time.
>
> > That is discrete time: t = [0,dt,2*dt,...(N-1)*dt];
>
> again, i really don't give a rat's ass about dt.  it's just a positive
> and real scaling factor (or a measure of time or something else in
> terms of an anthropocentric definition of the unit) and it doesn't
> change the phase angle.  dt could just as well be 1.

I wrote the exact formula because I am confused by your introduction
of the adjective "continuous". Now I am confused about the "dt" rant.

> > > however, i *do* worry about the inherent periodicity of the DFT.

Again, it's right there in the formula. Why the worry?

Greg