2
$\begingroup$

I am trying to practically determine filter group delay so that I can try to reduce or adjust for it. I have tried all kinds of methods but still can't get the right result. I know that theoretically to get the group delay you should differentiate the phase response with respect to frequency. How can I use this theoretical result to determine what delay will actually occur to my whole signal envelope? I ask because the theory does not seem to apply in practice in the way I have understood it.

I have tried to align filtered and unfiltered signal to see an actual shift but I do not seem to see any delay. I have tried adding a spike (Ok which is zero frequency) just so that I can compare the position of the spike in the original and the filtered signal, but the spike does not move in the filtered signal.

Additionally, some filters have both negative and positive theoretical group delay. In those filters should I take the average across all frequencies to determine the effective group delay?

I have added an example below in MATLAB. I used a linear phase FIR filter (a comb filter) which should have a delay of 20 samples according to MATLAB across all frequencies. But I have looked at the plots with the filtered and original signal on the same axis but could not get 20 samples at several frequencies.

%% Filter
b=[1,zeros(1,39),-1];%y(n)=x(n)-x(n-40)
a=1;

subplot(3,1,1)
grpdelay(b,1)

% Simulate
Fs=1000;
t=0:1/Fs:(5-1/Fs);
wi=blackman(length(t))';

spike=zeros(1,length(t));
spike(300)=0.02;%Feature
spike(150)=0.02;%feature

x1=sin(2*pi*49*t).*wi+spike;
x2=sin(2*pi*25*t).*wi;
x3=sin(2*pi*2*t).*wi;

x=[x1,x2,x3];

y=filter(cb,1,x);
subplot(3,1,2)
plot([x',y'])
title('x vs filtered x');
legend({'x','x-filtered'})

% show spike
subplot(3,1,3)
plot([x',y'])
title('x vs filtered x, zoomed to see spikes');
legend({'x','x-filtered'})
xlim([130,350])

The output of the above MATLAB code. I cannot determine where is the delay of 20 samples

My overall question is how can I practically measure the effective signal envelop delay due to filtering? I would like to match this measured delay with the theoretically calculated delay.

$\endgroup$
  • $\begingroup$ One method would be cross-correlation. Correlate the filter input and output, and look for the time lag that gives you the largest correlation value (in magnitude). That lag should correspond to the filter's group delay (20 samples in your case). $\endgroup$ – Jason R Jan 10 at 18:11
  • $\begingroup$ I was told off for doing cross correlation in a previous question because it is not general enough. dsp.stackexchange.com/questions/61964/… $\endgroup$ – Chika Jan 10 at 18:15
  • $\begingroup$ Yeah, that's a decent point: for an FIR filter, the correlation between the input and output at a particular time lag will be proportional to the magnitude of the corresponding filter coefficient. If you take the approach I suggested, then, you end up with a cross-correlation output that looks just like the magnitude of the filter taps. In that case, the group delay would be equal to the index of the largest tap. For symmetric (linear-phase) FIR filters, this is typically the case. $\endgroup$ – Jason R Jan 10 at 18:18
  • $\begingroup$ @Chika Did you try what I suggested at that other post; using a channel estimation approach since you have a complete replica of input and output? This should work best when your waveform has energy at all frequencies in your passband. You estimate the channel using a least squared approach and then from that you can get the group delay using the grpdelay command in Matlab, Octave or Python scipy.signal $\endgroup$ – Dan Boschen Jan 10 at 18:22
  • 1
    $\begingroup$ Does this answer your question? How to reliably compute the group delay of a comb filter $\endgroup$ – Dan Boschen Jan 12 at 16:38
4
$\begingroup$

The Least Mean Square solution to find the "channel" or response of the filter is provided by the following MATLAB/Octave Code using the input to the filter as tx and the output of the filter as rx. For more details on how this works, see this post: Compensating Loudspeaker frequency response in an audio signal:

function coeff = channel(tx,rx,ntaps)
    % Determines channel coefficients using the Wiener-Hopf equations (LMS Solution)
    % TX = Transmitted (channel input) waveform, row vector, length must be >> ntaps 
    % RX = Received (ch output) waveform, row vector, length must be >> ntaps 
    % NTAPS = Number of taps for channel coefficients
    % Dan Boschen 1/13/2020

    tx= tx(:)';   % force row vector
    rx= rx(:)';   % force row vector
    depth = min(length(rx),length(tx));
    A=convmtx(rx(1:depth).',ntaps);
    R=A'*A;       % autocorrelation matrix
    X=[tx(1:depth) zeros(1,ntaps-1)].';
    ro=A'*X;      % cross correlation vector
    coeff=(inv(R)*ro);   %solution
end

The case the OP uses of a comb filter is one of the most challenging as it depends on signal energy at each frequency for a solution (this is specifically why linear equalizers, which this function is doing if you swap rx and tx, do not perform well in frequency selective channels and result in noise enhancement at the null locations). Below the frequency response of the test filter as determined with MATLAB or Octave showing the multiple frequency nulls associated with such a comb filter:

b=[1,zeros(1,39),-1];
freqz(b,1,2^14)     % 2^14 samples to show nulls

Comb Filter Frequency Response

MATLAB or Octave script to demonstrate use of above function and determine the delay between output and input:

%% Filter with OP's example
b=[1,zeros(1,39),-1];     % numerator coefficients
a = 1;                    % denominator coefficients

%% Generate signal using OP's code
Fs=1000;
t=0:1/Fs:(5-1/Fs);
wi=blackman(length(t))';
rn=+rand(1,length(t))*.2;
x1=sin(2*pi*13*t).*wi +rn;
x2=sin(2*pi*25*t).*wi +rn;
x3=sin(2*pi*2*t).*wi +rn;
x=[x1,x2,x3];

% Filter
y=filter(b,a,x);

%% Test filter estimation
cf=channel(x,y,61);

%compare original and estimated channel
subplot(2,1,1)
stem(b)
title("Actual Channel Response")
xlabel("Sample Number")
subplot(2,1,2)

stem(cf)
title("Estimated Channel Response")
xlabel("Sample Number")

Actual and Estimated Channel Response

We could have used exactly 41 taps when calling the function channel and it would have resolved the solution properly; however it is best practice to start with a much longer filter length, evaluate the result and then reduce the taps accordingly. In actual practice under noise conditions using more taps than necessary will lead to noise enhancement, so it is good to do the final solution with the minimum taps necessary to capture the dominant tap weights.

Observe with the groupdelay plot using MATLAB and Octaves grpdelay command the issue with not being able to resolve the delay where no signal gets through the filter (it would be difficult for anything to determine the delay of a single tone at one of these frequencies that is nulled by the filter!) but it is able to accurately determined the delay where signal energy exists. Similarly the waveform itself must have energy at all frequencies where we are seeking a solution. The spectral density of the OP's test waveform was sufficiently spread out over all frequencies to be suitable for this purpose (this is the reason why pseudo-random waveform make for good "channel sounding" patterns).

This plot is for comparison to the OP's plot showing that the group delay for this filter is 20 samples.

Group Delay

Here is another test case showing how this works with a more reasonable channel (no deep nulls) with the following coefficients and frequency response below that:

b = [0.2 .4 -.3 .4 .3 .1];

![moderate channel distortion

The solution is indistinguishable between actual and estimate, so I added noise to make it more interesting, using same x and y from code above:

% add noise
noise = 0.351*randn(1,length(y));
yn = y + noise;
snr = 20*log10(std(yn)./std(noise));
%% Test filter estimation
cf=channel(x,yn,10);

Group Delay with SNR

$\endgroup$
  • 1
    $\begingroup$ For completeness please define a=1 before y=filter(b,a,x); in your code. $\endgroup$ – Chika Jan 15 at 12:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.