why fftshift() is useful.
in [DSP]
|
From: robert bristow-johnson on 27 Jun 2010 14:00 On Jun 27, 11:41 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > On Jun 26, 9:54 pm, robert bristow-johnson <r...(a)audioimagination.com> > wrote: > > > > > On Jun 24, 5:48 pm, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> > > wrote: > > > > >On Jun 24, 5:12=A0pm, "fisico32" <marcoscipioni1(a)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. > 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 = [0 1 2 3 4 5 6 7 8 7 6 5 4 3 2 1]'; X = fft(x); plot(angle(X + 1e-10)); versus this one: x = [8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7]'; X = 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
From: Greg Heath on 28 Jun 2010 01:12 On Jun 27, 2:00 pm, robert bristow-johnson <r...(a)audioimagination.com> wrote: > On Jun 27, 11:41 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > > On Jun 26, 9:54 pm, robert bristow-johnson <r...(a)audioimagination.com> > > wrote: > > > On Jun 24, 5:48 pm, "fisico32" <marcoscipioni1(a)n_o_s_p_a_m.gmail.com> > > > wrote: > > > > >On Jun 24, 5:12=A0pm, "fisico32" <marcoscipioni1(a)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
From: robert bristow-johnson on 28 Jun 2010 14:41 On Jun 28, 1:12 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > On Jun 27, 2:00 pm, robert bristow-johnson <r...(a)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 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. > > 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)); 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=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. > 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. however, i *do* worry about the inherent periodicity of the DFT. > 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. 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) >= 0 ==> angle(X) = 0 for the even indices (or odd indices in MATLABland). > real(X) < 0 ==> angle(X) = pi for the odd indices (or even inside of MATLAB). > However, what property of x guarantees that > > real(X) >= 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. 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. click the options for the thread. then hit "Fixed font". L8r, r b-j
From: Greg Heath on 28 Jun 2010 16:46 On Jun 28, 2:41 pm, robert bristow-johnson <r...(a)audioimagination.com> wrote: > On Jun 28, 1:12 am, Greg Heath <he...(a)alumni.brown.edu> wrote: > > On Jun 27, 2:00 pm, robert bristow-johnson <r...(a)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 > > what would that be? (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=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 > > 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. 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. and 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 = [0:dt:T] instead of t =[0:dt:T-dt] or tb = [-T/2:dt:T/2] instead of tb = [-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 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. > > > > 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)); > > 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=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. > > 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]; however, i *do* worry about the inherent > periodicity of the DFT. > > > 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. > > 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) >= 0 ==> angle(X) = 0 > > for the even indices (or odd indices in MATLABland). > > > real(X) < 0 ==> angle(X) = pi > > for the odd indices (or even inside of MATLAB). > > > However, what property of x guarantees that > > > real(X) >= 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. 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. > > click the options for the thread. then 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
From: Greg Heath on 28 Jun 2010 16:57 On Jun 28, 4:46 pm, Greg Heath <he...(a)alumni.brown.edu> wrote: > On Jun 28, 2:41 pm, robert bristow-johnson <r...(a)audioimagination.com> > wrote: -----SNIP > > > > 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. > > > click the options for the thread. then 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
|
Next
|
Last
Pages: 1 2 Prev: S'ware for Complex Test Signals Next: MATLAB indexing issue (was Re: Use of MATLAB fftshift) |