Measuring group delay through a linear phase filter

Question

I'm having a difficult time wrapping my head around a phase/ group delay problem and it corresponding to where you start sampling so I built up a system to test some of these issues. The system passes 2 sinusoids through a lowpass FIR filter with linear phase, then I take an NFFT length fft starting at some arbitrary point of my output data, then measure the phase of each sinusoid. The signal flow is as follows:

[Sinusoids] => [FIR (length 65)] => [FFT] => [phase measurement]

I'm trying to use the fact that group delay can be calculated from the phase information in the fft as:

but this is where I get lost because it seems that where you start sampling at the output would impact the phase measurement. i.e. you can't start at any arbitrary point. But in a black box system where you can't figure out what the correct output sample to start on is how can you accurately measure group delay? Below is my script. My group delay should be 32, but depending on the variable out_start_idx the group delay changes dramatically. Shouldn't it be that regardless of where you start sampling the output, the phase change between two tones should still be constant?

close all
clear all

fs = 1e6;

%GROUP DELAY of filter (N-1)/2 = 32
%since h = length 65
h = remez(2^6, [0 fs*.05 fs*.15 fs/2]/(fs/2), [1 1 0 0]); 

%arbitrarily large fft
NFFT = 2^19; 

bin_space  = fs/NFFT;

num_tones = 2; 
space = 1000; % spacing in bins
for nn = 1:num_tones
 freqs(nn) = space*(nn)*bin_space;
end

%random starting phases between -pi/2 pi/2
phi_start = rand(1,num_tones)*pi-pi/2; 
t = 0:999999; 

sig = 0; 
for nn = 1:num_tones
 sig = sig + exp(1i*2*pi*freqs(nn)*t/fs + phi_start(nn));
end

out = filter(h,1, [sig zeros(1,length(h)+ 10)]); 

%{
figure
subplot(211)
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(h,NFFT)))))
hold on
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(sig,NFFT)))))
hold off
subplot(212)
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(out,NFFT)))))
hold on
plot((-.5:1/NFFT:.5-1/NFFT)*fs, 20*log10(fftshift(abs(fft(h,NFFT)))))
hold off
%}
f = 0:fs/NFFT:fs-fs/NFFT;

% take NFFT amount of samples from arbitrary starting point 
% in output
out_start_idx = 1; 
data_test = out(out_start_idx:out_start_idx+NFFT-1); 

% get index of tones
idx_tones = (space:space:space*(num_tones)) + 1 

% take fft of data
fft_out = fft(data_test,NFFT);

% get complex output from each tone bin
tone_IQ = fft_out(idx_tones); 

% calculate phase of each tone
phi_tones = angle(tone_IQ)*180/pi;

% calculate group delay
phi_2 = diff(fliplr(phi_tones))
phi_1 = diff(fliplr(phi_start)*180/pi); 
tau = (-1/360) * (phi_2 - phi_1)./diff(freqs)
grp_delay = tau*fs

% maybe do some extra processing
phi_out = phi_tones; 

% linear fit calculated phases
p = polyfit(freqs,phi_out,1);
y = polyval(p, freqs); 

% calculate deviation from linear phase
deviation_from_linear_phase = phi_out -y; 

% plot deviation / etc
figure
plot(freqs, y)
hold on
plot(freqs, phi_out, 'r*')

figure
stem(deviation_from_linear_phase)

Answer 1

phi_start = rand(1,num_tones)*pi-pi/2;

The line above illustrates your problem. Phase measurements always need to referenced back to something. An absolute phase measurement without a well defined reference is meaningless.

I believe that the method outlined above that you don't measure the absolute phase but the phase of the output relative to the phase of the input. It probably assumes that you can observe both the input and the output. You can measure at an arbitrary time as long as you measure the input and the output simultaneously .

Answer 2

The reason the two tone method works to measure group delay as given is because group delay is the negative derivative of phase with respect to frequency. Since the system is linear phase (constant delay at all frequencies) then we can over any two frequencies determine the delay from the difference in phase between any two frequencies at the input and the the difference in phase between those two same frequencies at the output and from that determine the additional phase slope added due to the linear phase system between input and output.

The further in separation of the two frequencies, the more accurately we can estimate the slope in the presence of noise but otherwise given the linear slope this will work with any two frequencies chosen.