[Chapter 5.1] Meaning of General Linear Phase or Group Delay of a filter?

For me, one better description of the general linear phase (GLP) is a constant group delay. By definition, group delay is a negative derivative of phase [1] (In fact, the detailed phase/delay relationship can be derived but I suggest we just accept them.). Then, the Derivative (group delay) of a linear (phase) is constant and vise versa.

Then, let's move on to the meaning of group delay.
the group delay of a frequency represents the delay unit of the filter to that frequency.
So, the filter may treat different frequencies with different delay units. For an extremely bad example of a non-linear filter, the input signal 'do re mi' may becomes 're mi do' in the output. A GLP filter can guarantee that such weird condition will never happen.

Github link


-----------------------------------------------------------------------------------
Here I wrote an example. The example comes from Chapter 5.1.2 in the bible of DSP[2] and I just implemented the sample.

First, an IIR Filter with phase response like this is given.


Fig. 1.

Here is the group delay (a negative derivative of phase response) and the magnitude response. Please note that I delay the frequency in +-0.2 pi for about 150 units. BTW, the filter is a low pass filter, so a signal higher than 0.8pi is expected to be filtered.

Fig. 2.

Then, let's input a test signal like 'do re mi'. The signal x[n] are '0.8pi, 0.2pi, 0.4pi' in order. The corresponding frequency response is also provided.


Fig. 3.

And here is the output signal. The signal becomes 'empty, 0.4pi, 0.2pi'. The signal component with 0.8pi is filtered out as expected.



Fig. 4.

To make things more clear, here I point the number of Fig. 3 and Fig. 4 together. For the 0.4pi component, the group delay is about 6.39 units but the group delay of the 0.2pi component is about *153* unit. The output signal can confirm the prediction from the group delay response. That is why the 0.2pi component becomes that last in the output.



Fig. 5.


In summary,

1. linear phase equals constant group delay.

2. GLP FIR filter can guarantee such a scenario will never happen. But IIR can never achieve GLP. (But with the same mag frequency spectrum requirement, IIR usually can achieve the spec with lower delay (but not constant) in comparison with FIR.)

----------------------------------------------------------------------------------------------------------------------
Reference:


1. https://en.wikipedia.org/wiki/Group_delay_and_phase_delay

2. A. Oppenheim and R. Schafer, Discrete-Time Signal Processing 3rd. 2009

-----------------------------------------------------------------------------------
Matlab code

    


    %% System
    % H1[z]
    b1 = conv([1 -.98*exp(j*.8*pi)],[1 -.98*exp(-j*.8*pi)]);
    a1 = conv([1 -.8*exp(j*.4*pi)],[1 -.8*exp(-j*.4*pi)]);
    H1 = tf(b1,a1,-1,'Variable','z^-1');
 
    % H2[z]
    H2 = tf(1,1,-1,'Variable','z^-1');
    for k = 1:4
        ck = 0.95*exp(j*(0.15*pi+0.02*pi*k));
        ck_conj = conj(ck);
        b_tmp = conv([ck_conj -1],[ck -1]);
        b_tmp = conv(b_tmp,b_tmp);
        a_tmp = conv([1 -1*ck],[1 -1*ck_conj]);
        a_tmp = conv(a_tmp,a_tmp);
        H_tmp = tf(b_tmp,a_tmp,-1,'Variable','z^-1');
 
        H2 = series(H2,H_tmp);
    end

    % H[z]
    H = series(H1,H2);
 
    % Zero-Pole Plot, Fig. 5.2
    [b_h,a_h] = tfdata(H );
    b_h = cell2mat(b_h);
    a_h = cell2mat(a_h);

    figure;
    zplane(b_h,a_h);
    suptitle('Zero-Pole Plot, Fig 5.2');

    % System Response.
    L=1000;
    dw=2*pi/L;
    w = -pi:dw:pi-dw;
    HH=freqz(b_h,a_h,w);
    mag=abs(HH);
    phase=angle(HH);

    % Fig. 5.3
    figure;
    subplot(2,1,1);
    plot(w,phase);
    xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
    xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
    xlim([-pi pi]);
    yticks([-4 -2 -0 2 4]);
    ylabel('ARG[H(e^(^j^w^)]');
    xlabel('w');
    title('Phase response');

    subplot(2,1,2);
    plot(w,unwrap(phase));
    xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
    xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
    xlim([-pi pi]);
    ylabel('arg[H(e^(^j^w^)]');
    xlabel('w');
    title('Unwrap Phase response');
    suptitle('ARG/arg Plot, Fig 5.3');

    % Fig. 5.4
    figure;
    subplot(2,1,1);
    plot(w(1:end-1),-1*diff(unwrap(phase))./diff(w));
    xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
    xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
    xlim([-pi pi]);
    ylabel('grd[H(e^(^j^w^)]');
    title('Group Delay');

    subplot(2,1,2);
    plot(w,mag);
    xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
    xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
    xlim([-pi pi]);
    ylabel('|H(e^(^j^w^)|');
    title('Magnitude response');
    suptitle('GD/mag Plot, Fig 5.4');

    %% Signal
    M = 60;
    n = 0:M;
    w = 0.54-0.46*cos(2*pi*n/M);

    N = 512;
    x1 = zeros(1,N);
    x2 = zeros(1,N);
    x3 = zeros(1,N);
    dw = 2*pi/N;
    w_freq = -pi:dw:pi-dw;

    for i = 0:M
 
        x1(i+M) = w(i+1)*cos(0.2*pi*i);
        x2(i+2*M-1) = w(i+1)*cos(0.4*pi*i-pi/2);
        x3(i+1) = w(i+1)*cos(0.8*pi*i+pi/5);
    end
    x = x1+x2+x3;
    X = abs(fft(x));
    X = fftshift(X);

    % Fig. 5.5
    figure;
    subplot(2,1,1);
    plot(x);
    title('x[n]');
    xlim([0,300]);
    subplot(2,1,2);
    plot(w_freq,X);
    xticks([-pi -0.8*pi -0.6*pi -0.4*pi -0.2*pi 0 0.2*pi 0.4*pi 0.6*pi 0.8*pi pi]);
    xticklabels({'-\pi','-0.8\pi','-0.6\pi','-0.4\pi','-0.2\pi','0','0.2\pi','0.4\pi','0.6\pi','0.8\pi','\pi'});
    xlim([-pi pi]);
    ylabel('|H(e^(^j^w^)|');
    title('DTFT of X');
    suptitle('Input time/Freq., Fig 5.5');


    %% Output
    y = filter(b_h,a_h,x);

    % Fig. 5.6
    figure;
    plot(y);
    xlim([0,300]);
    xlabel('n');
    title('output y[n], Fig 5.6');


    %= Compre the Delay sample point.
    figure;
    subplot(2,1,1);
    plot(x);
    xlim([0,300]);
    xlabel('n');
    ylabel('x[n]');
    title('input');
    subplot(2,1,2);
    plot(y);
    xlim([0,300]);
    xlabel('n');
    ylabel('y[n]');
    title('output');

Comments

Popular posts from this blog

[Ch4 ADC/DAC] How to simulate ADC/DAC process in Matlab? How to actually reconstruct a signal nearly sampled in Nyquist Rate?

[Ch7. Filter Design] Why is it a bad idea to filter by zeroing out FFT bins? (為何不能DFT轉過去,直接設成想要的形狀,再IDFT轉回來?)