On Sat, 23 Nov 2013 12:25:57 +0100, Marcel M�ller <news.5.maazl@spamgourmet.org> wrote:>I want to analyze audio frequency response and group delay from room >measurements with microphone. Unfortunately I have no luck with the >group delay. Because of the cyclic nature of the phase I never know >whether neighborhood frequency bins have a positive or a negative delay. > >Current setup: >- Periodic reference signal with 2^18 samples length. >- Measure N periods on 2^18 samples and add the results. >- Calculate transfer function: complex division component by component > H[f] = FFT(microphone result[i]) / FFT(reference[i]) >- Calculate average over all delay by cross correlation (result is stable). >- Calculate delay per frequency by phase difference of neighborhood > frequencies > (arg(H[f]) - arg(H[f-1])) / delta f >- Find the delay value closest to the overall delay. > >This last one is the point that is often not unique. The complex room >response makes things complicated. Noise is no major problem. Repeated >measurements will show almost the same values. > >Any idea how to improve the result? > > >MarcelHi, Perhaps the following will be of some use to you: http://www.dsprelated.com/showarticle/69.php [-Rick-]
Calculate group delay from FFT data
Started by ●November 23, 2013
Reply by ●November 26, 20132013-11-26
Reply by ●November 26, 20132013-11-26
On 11/24/13 6:58 PM, Mac Decman wrote:> On Sat, 23 Nov 2013 19:39:13 -0800, robert bristow-johnson > <rbj@audioimagination.com> wrote: > >> >> >> >> alpha = 1.5/((2*pi)/3 - 4/pi) * 1/((frameSize/2)^2); % Peeters and Rodet would prefer 2.65 instead of 1.5 >> >> analWindow = linspace(-(frameSize-1)/2, (frameSize-1)/2, frameSize)'; >> analWindow = sqrt((frameSize/2)^2*alpha)*exp(-pi*alpha*(analWindow.^2)); % gaussian window >> >> ... >> xWindowed = x(round(n)+1:round(n)+frameSize) .* analWindow; % get and window the input audio frame >> >> X = fft(fftshift(xWindowed)); % do the FFT >> >> ... >> mag = log( abs(X(1:frameSize/2+1)) ); >> [maxmag, maxindex] = max(mag); >> delta = log( X(2:frameSize/2+1) ./ X(1:frameSize/2) ); >> realDelta = real(delta); >> realDeltaDelta = [0; diff(realDelta)]; % pad so that realDeltaDelta lines up with mag >> imagDelta = unwrap(imag(delta)); >> phaseDeltaDelta = unwrap(diff(imagDelta)); >> ... >> phaseDelta = cumsum([imagDelta(1); phaseDeltaDelta - (2*pi)*round(phaseDeltaDelta(maxindex-1)/(2*pi))]); >> phaseDelta = phaseDelta - (2*pi)*round(phaseDelta(maxindex))/(2*pi)); % principal value >> phase = cumsum([angle(X(1)); phaseDelta]); >> phase = phase - (2*pi)*round(phase(maxindex)/(2*pi)); % principal value > > hehehe this script is too dirty for me to run ;-)well, it appears that frameSize needs to be an odd integer, maybe not. looks sorta contradictory. and i dunno why i needed to round n. and i do not remember how i came up with the relationship between alpha and frameSize. this might be a teeny bit cleaner: alpha = 1.5/((2*pi)/3 - 4/pi) * 1/((frameSize/2)^2); % Peeters and Rodet would prefer 2.65 instead of 1.5 analWindow = linspace(-frameSize/2, frameSize/2 - 1, frameSize)'; analWindow = sqrt((frameSize/2)^2*alpha)*exp(-pi*alpha*(analWindow.^2)); % gaussian window ... xWindowed = x(n+1:n+frameSize) .* analWindow; % get and window the input audio frame X = fft( fftshift(xWindowed) ); % do the FFT ... mag = log( abs(X(1:frameSize/2+1)) ); [maxmag, maxindex] = max(mag); delta = log( X(2:frameSize/2+1) ./ X(1:frameSize/2) ); realDelta = real(delta); realDeltaDelta = [0; diff(realDelta)]; % pad so that realDeltaDelta lines up with mag imagDelta = unwrap(imag(delta)); phaseDeltaDelta = unwrap(diff(imagDelta)); ... phaseDeltaDelta = phaseDeltaDelta - (2*pi)*round(phaseDeltaDelta(maxindex-1)/(2*pi)); phaseDelta = cumsum([imagDelta(1); phaseDeltaDelta]); phaseDelta = phaseDelta - (2*pi)*round(phaseDelta(maxindex))/(2*pi)); % principal value phase = cumsum([angle(X(1)); phaseDelta]); phase = phase - (2*pi)*round(phase(maxindex)/(2*pi)); % principal value -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●November 26, 20132013-11-26
On Tue, 26 Nov 2013 11:48:24 -0800, robert bristow-johnson <rbj@audioimagination.com> wrote:>On 11/24/13 6:58 PM, Mac Decman wrote: >> On Sat, 23 Nov 2013 19:39:13 -0800, robert bristow-johnson >> <rbj@audioimagination.com> wrote: >> >>> >>> >>><snip> I think I may have tossed that one right past you. I always hate it when I'm working on a dsp program and I end up with lots of anal variables ;-) And you have the cumsum too ! Mark
Reply by ●November 26, 20132013-11-26
On 11/26/13 12:29 PM, Mac Decman wrote:> On Tue, 26 Nov 2013 11:48:24 -0800, robert bristow-johnson > <rbj@audioimagination.com> wrote: > >> On 11/24/13 6:58 PM, Mac Decman wrote: >>> On Sat, 23 Nov 2013 19:39:13 -0800, robert bristow-johnson >>> <rbj@audioimagination.com> wrote: >>> >>>> >>>> >>>> > <snip> > > I think I may have tossed that one right past you.maybe, but i saw the smiley face. problem is that i can't figger out from my own (old) code if frameSize was s'posed to be even or odd.> I always hate it > when I'm working on a dsp program and I end up with lots of anal > variables ;-) And you have the cumsum too ! >well, at least there isn't any [sp]rectal leakage. saw this in a psychiatrist's office: "You can't spell 'analysis' without 'anal'." glad he wasn't a gynecologist: https://www.youtube.com/watch?v=eQ87hBFrS_I if you haven't seen the movie, "slide to" the clip at 2:10 . -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●November 26, 20132013-11-26
On 26.11.13 15.23, Rick Lyons wrote:> Hi, > Perhaps the following will be of some use to you: > > http://www.dsprelated.com/showarticle/69.phpcool |;-) It is more convenient than dealing with the dirty phase unwrapping. Although I am unsure whether it is mathematically different in any way. Marcel
