DSPRelated.com
Forums

this fft code doesn't make sense

Started by BillJosephson November 26, 2006
BillJosephson skrev:
> % What do I fail to understand, when I expect this code to produce two > % waveforms of the same frequency, and identical fft peaks?
....
> % but the second one shows the peak at 121 Hz (which is 120).
No, it doesn't. The peak shows up in bin 121 which still is at 12 Hz. Look at the exponential term in your DFT. There you will find 2*pi*k*n/N These factors seem thrown together at random, but they are not: - The 2*pi factor is there because the DFT is defined in terms of angular frequencies. The n is a running index over samples in the input data sequence. - The N is the number of samples in the input data sequence. - The k is the bin index, and the factor you have trouble with understanding. To see how the k fits in the big bicture, we need first to draw the big picture. The DFT takes a signal that consist of N samples that are regularly sampled in time as input, and transform them to a signal that is regularly sampled in frequency. In time domain we have the sampling interval T, and the signal lasts for (N-1)T seconds. The signal duration changes as N changes. In frequency domain, we have the sampling frequency Fs = 1/T and the bin width dF = Fs/N. The spectrum is (N-1)*dF wide, or (N-1)/N*Fs. Here is the key point: The sampling frequency is fixed when N changes. The bin width dF changes when N changes. The question now is: At what frequency does coefficient X_k fall, given N and T? The k'th frequency occur at F_k = k*dF = k*Fs/N =k/(N*T). Looking at that expression, one might ask why the T factor is important or even necessary; and the answer is that it is neither. One can get rid of it by normalizing the frequency scale with respect to the sampling frequency: f_k = F_k/Fs = (k/N*T)/(1/T) = k/N And now we recognize the exonential term in the DFT: 2*pi*n*F_k -> 2*pi*n*k/N when the T factors are "distilled" away. To find out in what bin a certain frequency falls, you need to reverse this whole excercise: Given N and T, in what bin (indexed by k) contains the frequency f? First, find Fs and dF: Fs = 1/T dF = 1/(NT) Next, k is found by comparing ratios: f/Fs = k*dF/Fs f*T = (k/N*T)/(1/T) f*T = kT/NT f*T = k/N which, at last, gives k = f*T*N ======= and that's it. Just remember that k has to be an integer while the ratio f*T needs not be an integer. Rune
Ben Jackson wrote:
> On 2006-11-26, BillJosephson <billjosephson@hotmail.com> wrote: > > % now I have two fft results for a 12 hz sine wave. The first one > > correctly > > % shows the peak at 13 Hz (which is 12). > > % but the second one shows the peak at 121 Hz (which is 120). > > The first is a coincidence. Since your input sample was 1 second long, > the sample rate Fs and the size of the fft cancelled out. When you > increased the size of the sample and thus the fft they no longer cancel. > You'll see the same discrepancy in the vertical scale. Google fftplot() > (not sure if there's one in Matlab, there's not in Octave or Octave-forge > by default). It is a wrapper which should give you a plot in the correct > units (Hz vs power). > > -- > Ben Jackson AD7GD > <ben@ben.com> > http://www.ben.com/
Hi Ben, thanks for the message. Can you extricate your comment that a similar problem happens in the vertical scale? What do you mean? Thanks very much.
Rune Allnor wrote:
> BillJosephson skrev: > > % What do I fail to understand, when I expect this code to produce two > > % waveforms of the same frequency, and identical fft peaks? > .... > > % but the second one shows the peak at 121 Hz (which is 120). > > No, it doesn't. The peak shows up in bin 121 which still is at 12 Hz. > Look at the exponential term in your DFT. There you will find > > 2*pi*k*n/N > > These factors seem thrown together at random, but they are not: > > - The 2*pi factor is there because the DFT is defined in terms of > angular > frequencies. The n is a running index over samples in the input data > sequence. > - The N is the number of samples in the input data sequence. > - The k is the bin index, and the factor you have trouble with > understanding. > > To see how the k fits in the big bicture, we need first to draw the > big picture. > > The DFT takes a signal that consist of N samples that are > regularly sampled in time as input, and transform them to > a signal that is regularly sampled in frequency. In time domain > we have the sampling interval T, and the signal lasts for > (N-1)T seconds. The signal duration changes as N changes. > > In frequency domain, we have the sampling frequency Fs = 1/T > and the bin width dF = Fs/N. The spectrum is (N-1)*dF wide, > or (N-1)/N*Fs. > > Here is the key point: The sampling frequency is fixed when > N changes. The bin width dF changes when N changes. > > The question now is: At what frequency does coefficient X_k > fall, given N and T? > > The k'th frequency occur at F_k = k*dF = k*Fs/N =k/(N*T). > > Looking at that expression, one might ask why the T factor > is important or even necessary; and the answer is that it is > neither. One can get rid of it by normalizing the frequency > scale with respect to the sampling frequency: > > f_k = F_k/Fs = (k/N*T)/(1/T) = k/N > > And now we recognize the exonential term in the DFT: > > 2*pi*n*F_k -> 2*pi*n*k/N > > when the T factors are "distilled" away. > > To find out in what bin a certain frequency falls, you need > to reverse this whole excercise: > > Given N and T, in what bin (indexed by k) contains the > frequency f? > > First, find Fs and dF: > > Fs = 1/T > dF = 1/(NT) > > Next, k is found by comparing ratios: > > f/Fs = k*dF/Fs > f*T = (k/N*T)/(1/T) > f*T = kT/NT > f*T = k/N > > which, at last, gives > > k = f*T*N > ======= > > and that's it. Just remember that k has to be an integer > while the ratio f*T needs not be an integer. > > Rune
Rune, If you knew how much this message is helping me you'd know how much I appreciate your most kind help. I've been struggling with this off and on for the last year, and gotten help here and there, but this is the first message that is clear and complete enough that I think I have a chance of thorougly "getting it." Actually, I'm still working my way through it and implementing a function as I go. I may come back with a question, but at the moment am feeling like for once I have all the parts together with the plan all at the same time. Many, many thanks. Bill
dbell wrote:
> Bill, > > I think that what you don't understand is that one FFT is 10X longer > than the other but covers the same frequency range, so the index > spacing (in frequency) for the second case is 1/10 the index spacing of > the first case. > > So adjusting for the MATLAB index starting at one and not 0, you would > expect the second case peak to be at the same frequency as the first > case peak, which it is, but this corresponds to an index value that is > 10X the first case peak index, which it is. > > (121-1)=10*(13-1) > > Dirk Bell > DSP Consultant > > BillJosephson wrote: > > % What do I fail to understand, when I expect this code to produce two > > % waveforms of the same frequency, and identical fft peaks? > > > > t1 = (0:0.00048828:1); > > y1 = sin(2*pi*12*t1); > > t10 = (0:0.00048828:10); > > y10 = sin(2*pi*12*t10); > > > > % now I have two 12hz sine waves; y10 is 10 times as long as y1, and > > has 120 > > % oscillations instead of the 12 for y1.....otherwise the waves are > > % identical (I plotted and checked) > > > > f1=abs(fft(y1)); > > f10=abs(fft(y10)); > > > > % now I have two ffts, f10 is 10 times as long as f1 (and both the same > > as > > % y1 and y10) > > > > figure; plot(f1); figure; plot(f10); > > > > % now I have two fft results for a 12 hz sine wave. The first one > > correctly > > % shows the peak at 13 Hz (which is 12). > > % but the second one shows the peak at 121 Hz (which is 120). > > > > % thanks in advance for any further help (this is my second, and > > % hopefully better, question on this...)
dbell wrote:
> Bill, > > I think that what you don't understand is that one FFT is 10X longer > than the other but covers the same frequency range, so the index > spacing (in frequency) for the second case is 1/10 the index spacing of > the first case. > > So adjusting for the MATLAB index starting at one and not 0, you would > expect the second case peak to be at the same frequency as the first > case peak, which it is, but this corresponds to an index value that is > 10X the first case peak index, which it is. > > (121-1)=10*(13-1) > > Dirk Bell > DSP Consultant > > BillJosephson wrote: > > % What do I fail to understand, when I expect this code to produce two > > % waveforms of the same frequency, and identical fft peaks? > > > > t1 = (0:0.00048828:1); > > y1 = sin(2*pi*12*t1); > > t10 = (0:0.00048828:10); > > y10 = sin(2*pi*12*t10); > > > > % now I have two 12hz sine waves; y10 is 10 times as long as y1, and > > has 120 > > % oscillations instead of the 12 for y1.....otherwise the waves are > > % identical (I plotted and checked) > > > > f1=abs(fft(y1)); > > f10=abs(fft(y10)); > > > > % now I have two ffts, f10 is 10 times as long as f1 (and both the same > > as > > % y1 and y10) > > > > figure; plot(f1); figure; plot(f10); > > > > % now I have two fft results for a 12 hz sine wave. The first one > > correctly > > % shows the peak at 13 Hz (which is 12). > > % but the second one shows the peak at 121 Hz (which is 120). > > > > % thanks in advance for any further help (this is my second, and > > % hopefully better, question on this...)
Thank you, Dirk. I wonder if I would ever have figured out Matlab's odd (it is odd, not me, right?) convention about the first number being DC.
BillJosephson wrote:

   ...

> Thank you, Dirk. I wonder if I would ever have figured out Matlab's > odd (it is odd, not me, right?) convention about the first number being > DC.
The first bin is always DC (unless negative frequencies come first and DC is in the middle). What's odd -- unusual, that is,, is that it's numbered 1 instead of 0. None of the text material is written like that. Jerry -- Engineering is the art of making what you want from things you can get. &#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;
Rune Allnor wrote:

(snip)

> The DFT takes a signal that consist of N samples that are > regularly sampled in time as input, and transform them to > a signal that is regularly sampled in frequency. In time domain > we have the sampling interval T, and the signal lasts for > (N-1)T seconds. The signal duration changes as N changes.
I would have said NT. At least for periodic signals the period is NT.
> In frequency domain, we have the sampling frequency Fs = 1/T > and the bin width dF = Fs/N. The spectrum is (N-1)*dF wide, > or (N-1)/N*Fs.
This one I probably agree with, though it is hard to say that and still say that the FT is symmetric. The signal could be periodic in k-space. (snip) -- glen
BillJosephson skrev:
> Rune Allnor wrote: > > BillJosephson skrev: > > > % What do I fail to understand, when I expect this code to produce two > > > % waveforms of the same frequency, and identical fft peaks? > > .... > > > % but the second one shows the peak at 121 Hz (which is 120). > > > > No, it doesn't. The peak shows up in bin 121 which still is at 12 Hz. > > Look at the exponential term in your DFT. There you will find > > > > 2*pi*k*n/N > > > > These factors seem thrown together at random, but they are not: > > > > - The 2*pi factor is there because the DFT is defined in terms of > > angular > > frequencies. The n is a running index over samples in the input data > > sequence. > > - The N is the number of samples in the input data sequence. > > - The k is the bin index, and the factor you have trouble with > > understanding. > > > > To see how the k fits in the big bicture, we need first to draw the > > big picture. > > > > The DFT takes a signal that consist of N samples that are > > regularly sampled in time as input, and transform them to > > a signal that is regularly sampled in frequency. In time domain > > we have the sampling interval T, and the signal lasts for > > (N-1)T seconds. The signal duration changes as N changes. > > > > In frequency domain, we have the sampling frequency Fs = 1/T > > and the bin width dF = Fs/N. The spectrum is (N-1)*dF wide, > > or (N-1)/N*Fs. > > > > Here is the key point: The sampling frequency is fixed when > > N changes. The bin width dF changes when N changes. > > > > The question now is: At what frequency does coefficient X_k > > fall, given N and T? > > > > The k'th frequency occur at F_k = k*dF = k*Fs/N =k/(N*T). > > > > Looking at that expression, one might ask why the T factor > > is important or even necessary; and the answer is that it is > > neither.
******************************************************************** Forget about this for now! One can get rid of it by normalizing the frequency
> > scale with respect to the sampling frequency: > > > > f_k = F_k/Fs = (k/N*T)/(1/T) = k/N > > > > And now we recognize the exonential term in the DFT: > > > > 2*pi*n*F_k -> 2*pi*n*k/N > > > > when the T factors are "distilled" away.
*******************************************************************
> > To find out in what bin a certain frequency falls, you need > > to reverse this whole excercise: > > > > Given N and T, in what bin (indexed by k) contains the > > frequency f? > > > > First, find Fs and dF: > > > > Fs = 1/T > > dF = 1/(NT) > > > > Next, k is found by comparing ratios: > > > > f/Fs = k*dF/Fs > > f*T = (k/N*T)/(1/T) > > f*T = kT/NT > > f*T = k/N > > > > which, at last, gives > > > > k = f*T*N > > ======= > > > > and that's it. Just remember that k has to be an integer > > while the ratio f*T needs not be an integer. > > > > Rune > > > Rune, > > If you knew how much this message is helping me you'd know how much I > appreciate your most kind help. I've been struggling with this off and > on for the last year, and gotten help here and there, but this is the > first message that is clear and complete enough that I think I have a > chance of thorougly "getting it." Actually, I'm still working my way > through it and implementing a function as I go. I may come back with a > question, but at the moment am feeling like for once I have all the > parts together with the plan all at the same time. > > Many, many thanks.
Y're welcome. Seems like you would benefit from reading Rick Lyons' "Undertsanding Digital Signal Processing", Prentice Hall, 2004. After I posted yesterday I realized that one section does not fit in the overall post. I have outlined it above. The T factor needst to be preserved in the DFT all the way through, it can not be "distilled away" by mere arithmetics. The FFT is, however, a generic maths tool that is used in many applications where similar factors need not appear. Besides, leaving the T factor out also saves some precious FLOPs in computation-expensive applications. So by convention, everybody knows that the T factor ought to be present in the exponential term. I'm sorry for confusing you, but these conventions are second nature to me. I very seldom thing through the details about how things work. Rune
BillJosephson skrev:
> Rune Allnor wrote: > > BillJosephson skrev: > > > % What do I fail to understand, when I expect this code to produce two > > > % waveforms of the same frequency, and identical fft peaks? > > .... > > > % but the second one shows the peak at 121 Hz (which is 120). > > > > No, it doesn't. The peak shows up in bin 121 which still is at 12 Hz. > > Look at the exponential term in your DFT. There you will find > > > > 2*pi*k*n/N > > > > These factors seem thrown together at random, but they are not: > > > > - The 2*pi factor is there because the DFT is defined in terms of > > angular > > frequencies. The n is a running index over samples in the input data > > sequence. > > - The N is the number of samples in the input data sequence. > > - The k is the bin index, and the factor you have trouble with > > understanding. > > > > To see how the k fits in the big bicture, we need first to draw the > > big picture. > > > > The DFT takes a signal that consist of N samples that are > > regularly sampled in time as input, and transform them to > > a signal that is regularly sampled in frequency. In time domain > > we have the sampling interval T, and the signal lasts for > > (N-1)T seconds. The signal duration changes as N changes. > > > > In frequency domain, we have the sampling frequency Fs = 1/T > > and the bin width dF = Fs/N. The spectrum is (N-1)*dF wide, > > or (N-1)/N*Fs. > > > > Here is the key point: The sampling frequency is fixed when > > N changes. The bin width dF changes when N changes. > > > > The question now is: At what frequency does coefficient X_k > > fall, given N and T? > > > > The k'th frequency occur at F_k = k*dF = k*Fs/N =k/(N*T). > > > > Looking at that expression, one might ask why the T factor > > is important or even necessary; and the answer is that it is > > neither.
******************************************************************** Forget about this for now! One can get rid of it by normalizing the frequency
> > scale with respect to the sampling frequency: > > > > f_k = F_k/Fs = (k/N*T)/(1/T) = k/N > > > > And now we recognize the exonential term in the DFT: > > > > 2*pi*n*F_k -> 2*pi*n*k/N > > > > when the T factors are "distilled" away.
*******************************************************************
> > To find out in what bin a certain frequency falls, you need > > to reverse this whole excercise: > > > > Given N and T, in what bin (indexed by k) contains the > > frequency f? > > > > First, find Fs and dF: > > > > Fs = 1/T > > dF = 1/(NT) > > > > Next, k is found by comparing ratios: > > > > f/Fs = k*dF/Fs > > f*T = (k/N*T)/(1/T) > > f*T = kT/NT > > f*T = k/N > > > > which, at last, gives > > > > k = f*T*N > > ======= > > > > and that's it. Just remember that k has to be an integer > > while the ratio f*T needs not be an integer. > > > > Rune > > > Rune, > > If you knew how much this message is helping me you'd know how much I > appreciate your most kind help. I've been struggling with this off and > on for the last year, and gotten help here and there, but this is the > first message that is clear and complete enough that I think I have a > chance of thorougly "getting it." Actually, I'm still working my way > through it and implementing a function as I go. I may come back with a > question, but at the moment am feeling like for once I have all the > parts together with the plan all at the same time. > > Many, many thanks.
Y're welcome. Seems like you would benefit from reading Rick Lyons' "Undertsanding Digital Signal Processing", Prentice Hall, 2004. After I posted yesterday I realized that one section does not fit in the overall post. I have outlined it above. The T factor needst to be preserved in the DFT all the way through, it can not be "distilled away" by mere arithmetics. The FFT is, however, a generic maths tool that is used in many applications where similar factors need not appear. Besides, leaving the T factor out also saves some precious FLOPs in computation-expensive applications. So by convention, everybody knows that the T factor ought to be present in the exponential term. I'm sorry for confusing you, but these conventions are second nature to me. I very seldom thing through the details about how things work. Rune
BillJosephson skrev:
> Rune Allnor wrote: > > BillJosephson skrev: > > > % What do I fail to understand, when I expect this code to produce two > > > % waveforms of the same frequency, and identical fft peaks? > > .... > > > % but the second one shows the peak at 121 Hz (which is 120). > > > > No, it doesn't. The peak shows up in bin 121 which still is at 12 Hz. > > Look at the exponential term in your DFT. There you will find > > > > 2*pi*k*n/N > > > > These factors seem thrown together at random, but they are not: > > > > - The 2*pi factor is there because the DFT is defined in terms of > > angular > > frequencies. The n is a running index over samples in the input data > > sequence. > > - The N is the number of samples in the input data sequence. > > - The k is the bin index, and the factor you have trouble with > > understanding. > > > > To see how the k fits in the big bicture, we need first to draw the > > big picture. > > > > The DFT takes a signal that consist of N samples that are > > regularly sampled in time as input, and transform them to > > a signal that is regularly sampled in frequency. In time domain > > we have the sampling interval T, and the signal lasts for > > (N-1)T seconds. The signal duration changes as N changes. > > > > In frequency domain, we have the sampling frequency Fs = 1/T > > and the bin width dF = Fs/N. The spectrum is (N-1)*dF wide, > > or (N-1)/N*Fs. > > > > Here is the key point: The sampling frequency is fixed when > > N changes. The bin width dF changes when N changes. > > > > The question now is: At what frequency does coefficient X_k > > fall, given N and T? > > > > The k'th frequency occur at F_k = k*dF = k*Fs/N =k/(N*T). > > > > Looking at that expression, one might ask why the T factor > > is important or even necessary; and the answer is that it is > > neither.
******************************************************************** Forget about this for now! One can get rid of it by normalizing the frequency
> > scale with respect to the sampling frequency: > > > > f_k = F_k/Fs = (k/N*T)/(1/T) = k/N > > > > And now we recognize the exonential term in the DFT: > > > > 2*pi*n*F_k -> 2*pi*n*k/N > > > > when the T factors are "distilled" away.
*******************************************************************
> > To find out in what bin a certain frequency falls, you need > > to reverse this whole excercise: > > > > Given N and T, in what bin (indexed by k) contains the > > frequency f? > > > > First, find Fs and dF: > > > > Fs = 1/T > > dF = 1/(NT) > > > > Next, k is found by comparing ratios: > > > > f/Fs = k*dF/Fs > > f*T = (k/N*T)/(1/T) > > f*T = kT/NT > > f*T = k/N > > > > which, at last, gives > > > > k = f*T*N > > ======= > > > > and that's it. Just remember that k has to be an integer > > while the ratio f*T needs not be an integer. > > > > Rune > > > Rune, > > If you knew how much this message is helping me you'd know how much I > appreciate your most kind help. I've been struggling with this off and > on for the last year, and gotten help here and there, but this is the > first message that is clear and complete enough that I think I have a > chance of thorougly "getting it." Actually, I'm still working my way > through it and implementing a function as I go. I may come back with a > question, but at the moment am feeling like for once I have all the > parts together with the plan all at the same time. > > Many, many thanks.
Y're welcome. Seems like you would benefit from reading Rick Lyons' "Undertsanding Digital Signal Processing", Prentice Hall, 2004. After I posted yesterday I realized that one section does not fit in the overall post. I have outlined it above. The T factor needst to be preserved in the DFT all the way through, it can not be "distilled away" by mere arithmetics. The FFT is, however, a generic maths tool that is used in many applications where similar factors need not appear. Besides, leaving the T factor out also saves some precious FLOPs in computation-expensive applications. So by convention, everybody knows that the T factor ought to be present in the exponential term. I'm sorry for confusing you, but these conventions are second nature to me. I very seldom thing through the details about how things work. Rune