1
$\begingroup$

I'm trying to filter a 400 samples signal with various bandstop FIR filters (constant group delay), one at a time to see which one gives me the desired result. Every filter was build using the Kaiser Window method in Matlab. Each filter has a different number of coefficients. The smaller one has 782 and the biggest one 2460. I'm performing the filtering via fast convolution using the following Matlab code:

temp = ifft( fft(filterCoeffs,ffsSize)' .* fft( signalSamples, ffsSize) );

FilteredSignal = temp(offset+1:offset+length(signalSamples));

I'm trying to identify where inside the FilteredSignal vector of ffsSize are the 400 filtered samples. So far depending on the value of ffsSize and the amount of filter coefficients is the offset that I have to impose to recover the original 400 samples filtered. Even though I already know the right offset (where the useful filtered signal starts) I want to know the theoretical explanation of why there. I also want to know why is the group delay 1/2 of the amount of filter coefficients?

By the way if I use the Matlab function:

filter(filterCoeffs,1,signalSamples)

It returns 400 samples but not the right ones. Any idea why? Is like it doesn't identify the right offset

I have 2 cases:

No time aliasing (circular convolution equivalent to linear convolution)

ffsSize = length(filterCoeffs)+length(signalSamples) - 1;
aliasingSamples = 0
offset = groupDelay(Filter)+ length(signalSamples);

Without time aliasing

Time aliasing (circular convolution not equivalent to linear convolution)

ffsSize = length(filterCoeffs);
aliasingSamples = ffsSize - length(signalSamples);

if(groupDelay(Filter) < aliasingSamples) %this only happens when I use the 782 coefficient filter
    offset = groupDelay(Filter);
else
    offset = aliasingSamples;
end

with time aliasing

$\endgroup$
3
$\begingroup$

Lets put forward the intuition behind the concept of the group delay before further discussing how to find the delay of FIR filters.

Consider an input signal $x[n]$ of length $L_x$ which is nonzero between $n=0$ and $n=L_x-1$. And let a simplistic filter impulse response to be $h[n] = \delta[n-d]$. The output is immediately shown to be $y[n] = x[n-d]$ which is input shifted by $d$ samples to the right (for $d$ positive).

Let's look at the frequency response of the filter; $H(e^{j\omega}) = e^{-j \omega d}$ which is nothing but a phase term $\phi(\omega) = - \omega d$.

Consider a slightly complex filter $h[n] = 0.1 \delta[n] + 0.8\delta[n-1] + 0.1 \delta[n-2]$. Its frequency response can be written as $H(e^{j\omega}) = e^{-j 1 \omega} B(\omega)$ where $B(\omega) = 0.2 \cos(\omega n) + 0.8$ is a real valued function and has no phase-term. The output of this filter can be written in terms of delayed inputs due to each distinct impulse such as : $$y[n] = 0.1 x[n] + 0.8 x[n-1] + 0.1 x[n-2]$$

It can be proposed that $y[n]$ is dominated by the middle term $x[n-1]$ of a shift $1$ sample. (a more detailed analysis would say that low frequencies of $x[n]$ are dominated by the middle term and high frequencies are tending to cancel etc. which makes it a lowpass filter eventually)

In this latter example the phase function of the filter was $\phi(\omega) = -1\omega$ whereas the dominating component $x[n-1]$ at the output also had a delay of one sample to the left.

This observation can be used as an intuitive basis for the claim that the delay (or the group delay) associated with an FIR filter is given by the negative of the derivative of its phase response wrt frequency $\omega$; i.e., $$\boxed{ \tau = - \frac{d \phi(\omega)}{d\omega} }$$

Lets restrict the discussion to LTI, FIR symmetric filters which have linear phase responses. An LTI filter is described by its impulse response $h[n]$ which enables us to compute its response to an arbitrary input $x[n]$ via the discrete convolution sum:

$$ y[n] = x[n] \star h[n] = \sum_{k=-\infty}^{\infty} x[k] h[n-k] $$

where the dummy index $k$ would actually run over those valid range of nonzero multiplicands.

The filter also has a frequency response which is the DTFT of its impulse response: $$ H(e^{j\omega})= \mathcal{F} \{ h[n] \} =\sum_{n=-\infty}^{\infty} h[n]e^{-j\omega} = |H(e^{j\omega})| e^{j \phi_h(\omega)}$$

Now we define the group delay (in samples) of filter as: $$\tau = - \frac{d\phi_h(\omega)}{d\omega}$$

In general group delay will be a function of frequency $\omega$, however linear phase FIR filters have the property that their phase function is a linear function of frequency such as $$\phi_h(\omega) = K - d \omega$$

From which we deduce that the group delay will be $\tau = d$ which is independent of frequency; i.e. all components of an applied signal $x[n]$ will shift the same amount when filtered by a linear phase filter.

Given such a linear phase filter whose delay is $\tau = d$ and assume that you apply a very narrowband signal $x[n]$ such as: $$ x[n] = \cos(\omega_0 n) \frac{ \sin(\omega_1 n)}{\pi n} $$ where $\omega_1 \ll 1$ and $0 < \omega_0 < \pi$. Then the associated output for such a signal will be of the form $$ y[n]=A(\omega_0) x[n-d] = A(\omega_0) \cos(\omega_0 (n-d)) \frac{ \sin(\omega_1 (n-d))}{\pi (n-d)}$$

Finally, for example, it can be shown that a causal and symmetric FIR filter of length $M = 2L + 1$ from $n=0$ to $n=2L$ such that $h[n] = h[2L-n]$ will have a linear-phase term of the form $\phi_h(\omega) = K - \omega L$. Where $K$ can be zero. Hence its group delay will be $L$ samples. Lets show it here:

$$ \begin{align} H(e^{j\omega}) &= \sum_{n=0}^{2L} h[n] e^{-j \omega n} \\ &= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{n=L+1}^{2L} h[n] e^{-j \omega n} + h[L] e^{-j \omega L} \\ &= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{n=L+1}^{2L} h[2L-n] e^{-j \omega n} + h[L] e^{-j \omega L} \\ \end{align} $$ Substitute $2L-n = m$ in the second sum: $$ \begin{align} H(e^{j\omega}) &= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{m=L-1}^{0} h[m] e^{-j \omega (2L-m)} + h[L] e^{-j \omega L} \\ &= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + e^{-j\omega 2L} \sum_{m=L-1}^{0} h[m] e^{j \omega m} + h[L] e^{-j \omega L} \\ &= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + e^{-j\omega 2L} \sum_{m=0}^{L-1} h[m] e^{j \omega m} + h[L] e^{-j \omega L} \\ \end{align} $$ Recognize the similarity of the first sum and second sum such that denoting the first sum as $A(\omega)$ then the second sum is $A(\omega)^*$ (assuming real $h[n]$), Hence

$$ H(e^{j\omega}) = A(\omega) + e^{-j\omega 2L} A(\omega)^* + h[L] e^{-j \omega L} \\ $$

Taking all of them into the $e^{-j\omega L}$ parenthesis yields: $$ H(e^{j\omega}) = e^{-j \omega L} \left( A(\omega)e^{j \omega L} + A(\omega)^* e^{-j \omega L} + h[L] \right) \\ $$

Since $A(\omega)e^{j \omega L} + A(\omega)^* e^{-j \omega L} = 2 \mathcal{Re} \{ A(\omega)e^{j \omega L} \}$

We reduce the frequency response of the FIR filter to: $$H(e^{j\omega}) = e^{-j \omega L} ( 2 \mathcal{Re} \{ A(\omega)e^{j \omega L} \} + h[L] ) $$

Where the thing in the parenthesis is a real valued function and has zero phase term, therefore we call it $B(\omega)$ and deduce that the DTFT $H(e^{j\omega})$ of the causal, symmetric, real valued impulse response $h[n]$ is

$$ \boxed{ H(e^{j\omega}) = B(\omega) e^{-j\omega L} }$$

from which we find that the associated phase response of the filter is: $$\phi_h(\omega) = -\omega L $$

and the associated group delay is therefore $$\boxed{ \tau = - \frac{d \phi_h(\omega)} {d\omega} = L}$$

samples...

% This script demonstrates the use of time-aliasing in an advantageous way
% in obtaining the "central" portion of the circular convolution implemented
% via the DFT approach...
%
% Let x[n] of length Lx and h[n] of length Lh = 2*L+1, are convolved linearly
% 
%                     y[n] = x[n] * h[n]
%  
% so that "full" length of y[n] is:  Ly = Lx + Lh -1
%
% However in some cases we are only interested in those samples of length Lx 
% that relate to input x[n] which are at the "central" portion of y[n], whereas
% the initial and final L samples carry transient information. Hence we take
% 
%                    yc[n] = y[L:L+Lx-1]    
%
% When a linear convolution is implemented with DFT, the length of DFT must
% satisfy the following; if we are to avoid time-aliasing: M >= Lz=Lx+Lh-1
% However this method computes ALL the full samples of y[n] in addition to the
% central ones. Hence if we wish to compute only the central portion of
% linear convolution y[n], then we can allow some aliasing in the DFT method
% The amount of which is indicated in the length of DFT as:
%
%                    M = Lz - L
%
% Then this will produce a result ya[n] whose last Lx samples carry the 
% required information such that 
%
%                    yc[n] = y[L:L+Lx-1] = ya[Lz-L-Lx:Lz-L-1]
%
% in matlab notation the indices will be added by 1 as:
%
%                    yc(:) = y(L+1:L+Lx) = ya(Lz-L-Lx+1:Lx-L);
%
% ==============================================================================
clc; clear all; close all;

% S0 - Define the parameters:
% ---------------------------
Lx = 32;
L = 15;
Lh = 2*L+1;
Lz = Lx+Lh-1;

% S1 - Generate the signals:
% --------------------------
x = randn(1,Lx);
h = randn(1,Lh);


% S2 - implement time-domain convolution first:
% ---------------------------------------------
yt = conv(x,h);
y1 = yt(L+1:L+Lx);           % get the central portion of the convolution


% S3 - Implement frequency domain DFT based aliased-circular conv:
% ----------------------------------------------------------------
M = Lz-L;       % DFT length (with aliasing)

yt = real( ifft( fft(x,M).*fft(h,M) , M) );     % compute aliased result
y2 = yt(L+1:end);

% S4 - compare the results
% ------------------------
figure,stem(y1,'b')
hold on
stem(y2,'g')

figure,stem(y1-y2);title('the error');
$\endgroup$
  • $\begingroup$ +1Thank you very much for your answer! However, I'm also interested in knowing the theory behind it. Would you mind also answering: Why is the group delay 1/2 of the amount of the FIR filter coefficients? Why is Matlab filter function not working properly? Why is the useful filtered signal at groupDelay(Filter)+ length(signalSamples) for linear convolution equivalent and groupDelay(Filter) or aliasingSamples for circluar convolution with aliasing? $\endgroup$ – VMMF Oct 26 '17 at 14:16
  • $\begingroup$ Ok let me add some theory, but not much.. $\endgroup$ – Fat32 Oct 26 '17 at 14:23

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.