Calculate group delay from FFT data

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?


Marcel

Unwrap  H(f ) first before taking phase derivative. If the unwrapped phase has large bin- to- bin  jumps that indicate that the unwrap failed, then you need a longer sequence (room impulse response may have a longer duration than your sequence length). 
You can also obtain the impulse response directly and then apply a window before the fft. 

Bob

Am 23.11.2013 12:25, schrieb Marcel M�ller:
> 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.

What's the spectral content of that? Is it mostly white?

> - 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])

Why don't you simply use P. Welch's method for transfer function
estimation? In Octave, you would type

  [transfer,frequency] = pwelch(reference,recorded,...
      window,overlap,Nfft,Fs,'trans','whole');

You could use a sine sweep or white noise for the reference signal. Make
sure that 'reference' and 'recorded' are somewhat aligned (the more the
better), otherwise transfer function estimation will be off.

The group delay is simply the derivative of the phase response. For an
appropriate large Nfft, you can estimate the derivative using a finite
difference. Don't forget to scale this according to the sampling rate
and Nfft.

After zeroing the coefficients in 'transfer' that could not be
determined because the frequencies are not part of the reference signal,
you could also directly invoke the ifft function to get an idea of the
impulse response.

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

I know. Lots of reflections here and there. Basically, it's a multipath
channel. I recently determined my room's impulse response (including
speaker+mic) by playing white noise + using pwelch --- just to gen an
idea of how reverberant it is.

Cheers!
SG

On 11/23/13 5:22 AM, radams2000@gmail.com wrote:
> Unwrap  H(f ) first before taking phase derivative.

you can also get derivatives of phase directly from the real and imag 
parts *before* calculating the phase with any arg() or arctan() 
function.  you can integrate that derivative to get the unwrapped phase. 
  a line of simple MATLAB code i have used is this:


numSamples = length(x);

X = fft(x);

phase = cumsum([angle(X(1)); angle( X(2:numSamples/2) ./ 
X(1:numSamples/2-1) ) ]);

but it still assumes that the change in phase less than +/- 180 degrees 
between neighboring bins.

i have experimented with modeling the phase as quadratically increasing 
with frequency (like what you would get with a linear "sweep" or chirp) 
and it looked something like this

diff_diff_phase = unwrap(diff(unwrap(( angle( X(2:numSamples/2) ./ 
X(1:numSamples/2-1) ))));
diff_phase = cumsum([angle(X(2)./X(1)); diff_diff_phase]);
phase = cumsum([angle(X(1));  diff_phase]);

i think that might work, but i could readily find that old code.  i'm 
still looking.




> If the unwrapped phase has large bin-to-bin jumps that indicate that the unwrap failed, then you need a longer sequence (room impulse response may have a longer duration than your sequence length).
> You can also obtain the impulse response directly and then apply a window before the fft.


i kinda would assume doing that in a spectrum analyzer that might 
display group delay.

i remember doing this in the past and being a little dissatisfied with 
the unwrap() function.

i found the code i was looking for.  here are snippets:

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)+phaseDelta(maxindex))/(4*pi));	% 
principal value
	phase = cumsum([angle(X(1)); phaseDelta]);	
	phase = phase - (2*pi)*round(phase(maxindex)/(2*pi));	% principal value



you might have to unwrap some lines of code here, but this is how i 
dealt with the phase unwrapping problem in MATLAB.

but group delay (using finite difference instead of derivative) can be 
calculated directly from the real and imag parts without needing to deal 
with phase unwrapping.  but i agree with Bob that you gotta unwrap phase 
*somehow* before applying the derivative to the phase to avoid some 
nasty spikes in the derivative.  i just think skipping phase and 
calculating its derivative directly from real and imag parts obviates 
the phase unwrapping problem for the same of group delay.  maybe you 
gotta unwrap a higher derivative, i just don't remember.

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

in copying (and modifying to exclude some provincial stuff) i messed up 
one line sorta obviously.



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




-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

Hi,

in a multipath environment, group delay (as a single number) is not defined
unambiguously. Results will depend on the algorithm.


Here is a code snippet that determines "the" delay by locating the
crosscorrelation peak between input and output:
http://www.dsprelated.com/showcode/288.php

It achieves subsample resolution by oversampling the signal.
If that's what you want, it's fairly robust, if computationally expensive
and not really smart as an algorithm.


I've got an alternative, phase-based algorithm, here:
http://www.dsprelated.com/showarticle/26.php
It's mainly used for radio frequency amplifiers, not recommended for audio
paths with highly varying (or ill-defined) group delay. 
Have a look at the blog for comments on phase unwrapping: Generally, it
fails catastrophically if the channel response has zeros or drops below the
noise floor, a single "bad" FFT bin is enough.	 

_____________________________		
Posted through www.DSPRelated.com

If you want the group delay _vs frequency_,  as on a network analyzer
instrument - it is straightforward to detect discontinuities in GD(f) and
correct the phase to minimize them.
Problem is, it doesn't fail gracefully if the signal gets too noisy. 

	 

_____________________________		
Posted through www.DSPRelated.com

On 23.11.13 18.59, sg wrote:
> Am 23.11.2013 12:25, schrieb Marcel M�ller:
>> 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.
>
> What's the spectral content of that? Is it mostly white?

Mostly yes. The weight function is f^-0,25, something between white and 
pink noise.

>> - 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])
>
> Why don't you simply use P. Welch's method for transfer function
> estimation? In Octave, you would type

Never heard of that. As far as I see it is mainly a method to estimate 
the transfer function from an (theoretically) infinite measurement with 
a infinite long reference (non-deterministic, e.g. a noise generator).

I have no need to go this long way round. My reference is deterministic 
and contains only energy at frequencies that do very exactly fit to the 
FFT bins. There is almost no energy at other frequencies. In fact it 
repeats after the FFT's length. So I get a very good SNR even when 
environment noise is there.


> The group delay is simply the derivative of the phase response. For an
> appropriate large Nfft, you can estimate the derivative using a finite
> difference.

I know. I hoped that 5.5 seconds FFT length were enough for this 
purpose. Maybe I was wrong. It turns out that there were some bugs in 
the code that reduced the maximum length. I fixed that and retried with 
2^19 FFT length and now the delay seems to be reasonably smooth. Only at 
very low frequencies (below 20 Hz) the response is still jumpy. But I do 
not care about that. Normally this frequencies are turned off.

> After zeroing the coefficients in 'transfer' that could not be
> determined because the frequencies are not part of the reference signal,
> you could also directly invoke the ifft function to get an idea of the
> impulse response.

I find it quite difficult to see any useful information from an impulse 
response. But I will check that.

>> - 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.
>
> I know. Lots of reflections here and there. Basically, it's a multipath
> channel.

Ack.


Marcel

On Sat, 23 Nov 2013 19:39:13 -0800, robert bristow-johnson
<rbj@audioimagination.com> wrote:

>
>
>
>
>in copying (and modifying to exclude some provincial stuff) i messed up 
>one line sorta obviously.
>
>
>
>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 ;-)

Mark

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

I'm not sure if I was looking for a meaningful group delay number
(what ever you want that to mean for analysing a room) that I would be
setting up the measurement system as you have.

I think that I would first use some type of deterministic stimulus
which has an IF curve weighted to excite the room modes properly.  I
would then use some type of decomposition method on the recovered data
to break the response into analytic signal components.


Mark.