DSPRelated.com
Forums

FFTs of FFTs

Started by jaso...@gmail.com June 24, 2006
>Hi jason, > >I think ur friend may be referring to the technique published in 2001 >DAFx Conference. Please note this technique is not the same as >cepstrum. Both cepstrum and 'Fourier of Fourier Tx' can however be >related to pitch >You may want to look at this paper if this is what you are actually >referring to >http://www.csis.ul.ie/dafx01/proceedings/papers/marchand.pdf > >Regards, >Raghuram
ok, I tried to implement the algorithm in that paper, but somehow it doesn't seem to work ... I'm still unsure about the whole fft-shifting stuff ... can somebody maybe have a look at the code and tell me where I'm wrong?
>> S = flipud(specgram(signal, 512, 11025, hamming(512))); >> S = abs((S .* conj(S) / 512)).^2; >> >> temp = []; >> for i=1:size(S,2) >> curr = fftshift(S(:,i)); >> curr = curr - mean(curr); >> curr = curr(128:257); >> S2 = fftshift(fft(curr, 2048)); >> S2 = abs(S2 .* conj(S2) / 512) .^ 2; >> S2 = S2(1025:2048); >> >> temp = [temp min(find(S2 == max(S2)))]; >> end
MrBurnst wrote:

> can somebody maybe have a look at the code and tell me where I'm > wrong? > >>> S = flipud(specgram(signal, 512, 11025, hamming(512))); >>> S = abs((S .* conj(S) / 512)).^2;
I don't know if this solves the problem, but S conj(S) = |S|^2 to numerical precision. I'd choose the latter way so the imaginary part comes out as stricly zero. Martin -- Quidquid latine scriptum sit, altum viditur.
>MrBurnst wrote: > >> can somebody maybe have a look at the code and tell me where I'm >> wrong? >> >>>> S = flipud(specgram(signal, 512, 11025, hamming(512))); >>>> S = abs((S .* conj(S) / 512)).^2; > >I don't know if this solves the problem, but S conj(S) = |S|^2 to >numerical precision. I'd choose the latter way so the imaginary part >comes out as stricly zero. > > >Martin > >-- >Quidquid latine scriptum sit, altum viditur. >
thanks, but I don't think that's the problem. As far as I can tell, the imaginary part is 0.... The algo just allways pics the first bin, since it is (by far) the highest ... I think I'm having trouble with the whole fft-shifting / dc-component stuff, which I haven't quite figured out ... anybody help?
Martin Eisenberg wrote:
> MrBurnst wrote: > > >>can somebody maybe have a look at the code and tell me where I'm >>wrong? >> >> >>>>S = flipud(specgram(signal, 512, 11025, hamming(512))); >>>>S = abs((S .* conj(S) / 512)).^2; > > > I don't know if this solves the problem, but S conj(S) = |S|^2 to > numerical precision. I'd choose the latter way so the imaginary part > comes out as stricly zero. > > > Martin >
I haven't been following this thread. *BUT!!!* " ... to numerical precision" *??????* if S = a + bi then | S|^2 = a^2 + b^2 and conj(S) = a - bi _and_ a + bi X a - bi --------------- -abi + b^2 a^2 + abi ------------- a^2 + 0i + b^2 So where is " ... to numerical precision" *?* ;) So I never got my BSEE, but ... ;} ;] ;) OWL _ducks_ for cover
Richard Owlett wrote:
> Martin Eisenberg wrote: >> MrBurnst wrote: >> >> >>> can somebody maybe have a look at the code and tell me where I'm >>> wrong? >>> >>>>> S = flipud(specgram(signal, 512, 11025, hamming(512))); >>>>> S = abs((S .* conj(S) / 512)).^2; >> >> >> I don't know if this solves the problem, but S conj(S) = |S|^2 to >> numerical precision. I'd choose the latter way so the imaginary part >> comes out as stricly zero. >> >> >> Martin >> > > I haven't been following this thread. > > *BUT!!!* > > " ... to numerical precision" *??????* > > if S = a + bi then | S|^2 = a^2 + b^2 > and conj(S) = a - bi > _and_ > > a + bi > X a - bi > --------------- > -abi + b^2 > a^2 + abi > ------------- > a^2 + 0i + b^2 > > So where is " ... to numerical precision" *?* ;) > So I never got my BSEE, but ... ;} ;] ;) > > OWL _ducks_ for cover
People ask all the time why a filter coefficient that theory supposes to be zero comes out to 3.084e-19, and wonder what's wrong with the routine that calculates it. Would you bet much that the values of -iab and +iab, separately calculated with floating point will sum to precisely zero? No? I thought not. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Jerry Avins wrote:

> Richard Owlett wrote: > >> Martin wrote: >> >>> MrBurnst wrote: >>> >>> >>>> can somebody maybe have a look at the code and tell me where I'm >>>> wrong? >>>> >>>>>> S = flipud(specgram(signal, 512, 11025, hamming(512))); >>>>>> S = abs((S .* conj(S) / 512)).^2; >>> >>> >>> >>> I don't know if this solves the problem, but S conj(S) = |S|^2 to >>> numerical precision. I'd choose the latter way so the imaginary part >>> comes out as stricly zero. >>> >>> >>> Martin >>> >> >> I haven't been following this thread. >> >> *BUT!!!* >> >> " ... to numerical precision" *??????* >> >> if S = a + bi then | S|^2 = a^2 + b^2 >> and conj(S) = a - bi >> _and_ >> >> a + bi >> X a - bi >> --------------- >> -abi + b^2 >> a^2 + abi >> ------------- >> a^2 + 0i + b^2 >> >> So where is " ... to numerical precision" *?* ;) >> So I never got my BSEE, but ... ;} ;] ;) >> >> OWL _ducks_ for cover > > > People ask all the time why a filter coefficient that theory supposes to > be zero comes out to 3.084e-19, and wonder what's wrong with the routine > that calculates it. Would you bet much that the values of -iab and +iab, > separately calculated with floating point will sum to precisely zero? > No? I thought not. > > Jerry
*OH NO* I'm agreeing with Jerry ;} I think you state the obverse to my perverse. Mr. Eisenberg should have said S conj(S) = |S|^2 which *may be* approached to limits of "numerical precision".
>> >> >> People ask all the time why a filter coefficient that theory supposes
to
>> be zero comes out to 3.084e-19, and wonder what's wrong with the
routine
>> that calculates it. Would you bet much that the values of -iab and
+iab,
>> separately calculated with floating point will sum to precisely zero? >> No? I thought not. >> >> Jerry > >*OH NO* I'm agreeing with Jerry ;} > >I think you state the obverse to my perverse. > >Mr. Eisenberg >should have said S conj(S) = |S|^2 >which *may be* approached to limits of "numerical precision". > >
Somehow I don't believe that differences in that order of magnitude (e-19) are responsable for the behavior I'm observing. I have a really strong (most of the time THE strongest in the spectrum) element in the 0est bin. What I'm doing is taking the FFT of a FFT. Now, the first one has 512 bins, and the resulting spectrum is zero-padded to a length of 2048. Can this procedure somehow introduce some kind of offset? As I said, first, I just use matlab's specgram method, which outputs one half of the complex spectrum. On each window of the spectrum, I comute the magnitude spectrum and apply the FFT, which in turn returns a complex spectrum, of which I again compute the power spectrum, but only using the first 1024 bins.... (??)
MrBurnst wrote:
>>> >>> People ask all the time why a filter coefficient that theory supposes > to >>> be zero comes out to 3.084e-19, and wonder what's wrong with the > routine >>> that calculates it. Would you bet much that the values of -iab and > +iab, >>> separately calculated with floating point will sum to precisely zero? >>> No? I thought not. >>> >>> Jerry >> *OH NO* I'm agreeing with Jerry ;} >> >> I think you state the obverse to my perverse. >> >> Mr. Eisenberg >> should have said S conj(S) = |S|^2 >> which *may be* approached to limits of "numerical precision". >> >> > > Somehow I don't believe that differences in that order of magnitude (e-19) > are responsable for the behavior I'm observing. I have a really strong > (most of the time THE strongest in the spectrum) element in the 0est bin. > What I'm doing is taking the FFT of a FFT. Now, the first one has 512 > bins, and the resulting spectrum is zero-padded to a length of 2048. > Can this procedure somehow introduce some kind of offset? As I said, > first, I just use matlab's specgram method, which outputs one half of the > complex spectrum. On each window of the spectrum, I comute the magnitude > spectrum and apply the FFT, which in turn returns a complex spectrum, of > which I again compute the power spectrum, but only using the first 1024 > bins.... (??)
If you tack all the zeros at the end, your real samples are in the first quarter of the window. That's bound to bias _something_. (Hardly any operation is trivial, even id powerful tools make it seem so.) Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
>> >> Somehow I don't believe that differences in that order of magnitude
(e-19)
>> are responsable for the behavior I'm observing. I have a really strong >> (most of the time THE strongest in the spectrum) element in the 0est
bin.
>> What I'm doing is taking the FFT of a FFT. Now, the first one has 512 >> bins, and the resulting spectrum is zero-padded to a length of 2048. >> Can this procedure somehow introduce some kind of offset? As I said, >> first, I just use matlab's specgram method, which outputs one half of
the
>> complex spectrum. On each window of the spectrum, I comute the
magnitude
>> spectrum and apply the FFT, which in turn returns a complex spectrum,
of
>> which I again compute the power spectrum, but only using the first
1024
>> bins.... (??) > >If you tack all the zeros at the end, your real samples are in the first
>quarter of the window. That's bound to bias _something_. (Hardly any >operation is trivial, even id powerful tools make it seem so.) > >Jerry >-- >Engineering is the art of making what you want from things you can get. >����������������������������������������������������������������������� >
Well, that's according to the procedure proposed in http://www.csis.ul.ie/dafx01/proceedings/papers/marchand.pdf changing those numbers (the padding) has no positive effects...
MrBurnst wrote:
>>> >>>People ask all the time why a filter coefficient that theory supposes > > to > >>>be zero comes out to 3.084e-19, and wonder what's wrong with the > > routine > >>>that calculates it. Would you bet much that the values of -iab and > > +iab, > >>>separately calculated with floating point will sum to precisely zero? >>>No? I thought not. >>> >>>Jerry >> >>*OH NO* I'm agreeing with Jerry ;} >> >>I think you state the obverse to my perverse. >> >>Mr. Eisenberg >>should have said S conj(S) = |S|^2 >>which *may be* approached to limits of "numerical precision". >> >> > > > Somehow I don't believe that differences in that order of magnitude (e-19) > are responsable for the behavior I'm observing. I have a really strong > (most of the time THE strongest in the spectrum) element in the 0est bin. > What I'm doing is taking the FFT of a FFT. Now, the first one has 512 > bins, and the resulting spectrum is zero-padded to a length of 2048. > Can this procedure somehow introduce some kind of offset? As I said, > first, I just use matlab's specgram method, which outputs one half of the > complex spectrum. On each window of the spectrum, I comute the magnitude > spectrum and apply the FFT, which in turn returns a complex spectrum, of > which I again compute the power spectrum, but only using the first 1024 > bins.... (??)
*I AGREE!* ; I was pedantically pointing out Mr. Eisenberg's analysis faux pas.