Huo Jiaquan wrote:

> I am reading a piece of C code supposedly implementing a 1st
> order all pass filter of transfer function
> A(z) = (C+z^-1)/(1+C*z^-1) with the output decimated by a factor
> of 2

More precisely, it applies A(z) to the even-indexed input subsequence 
and writes out the result with stride two. Let me restate the loop 
more succinctly:

for(i = 0;  i < N/4;  ++i) {
  t = in[4*i+0] - C*d;
  out[4*i+0] = d + C*t;
  d = in[4*i+2] - C*t;
  out[4*i+2] = t + C*d;
}

This is equivalent to the loop:

for(i = 0;  i < N/2;  ++i) {
  t = in[2*i] - C*d;
  out[2*i] = d + C*t;
  d = t;
}

> I made the following derivation:
>      Y(z) = A(z)X(z).
> Since we want to do decimation, I multiply the numerator and
> denumerator by (1-Cz^-1). Note that the conjugate on C is
> ignored. We get 
>      Y(z) = (1 - Cz^-1)(C + z^-1) / (1 - C^2z^-2)X(z).
> Moving the denumerator to the left and then moving z^-2Y(z) term
> to the right, it follows
>      Y(z) = CX(z) + (1 - C^2)z^-1X(z) - Cz^-2X(z) - C^2z^-2Y(z).

The last term above and in Q(z) below should have positive sign.

> To avoid buffering Y(z) (as you can see, it's named temp1 in the
> code, it's just an intermidiate result), I group all terms with
> delay together, call it Q(z)
>      Q(z) = (1 - C^2)z^-1X(z) - Cz^-2X(z) - C^2z^-2Y(z).
> Let
>      I(z) = Q(z) / (1 - C^2),
> we have
>      Y(z) = I(z)+C[X(z)-CI(z)],
> which is the first 2 lines of the code, if data0=I(z).

Right.

> But then the 3rd and 4th lines are very odd.

Let's ignore them for a second and carry on with your analysis.

> For updating I(z) from z^-2I(z), substitute Y(z) as in the above
> equation into the definition of I(z), one can get
>       I(z) = z^-1X(z) - Cz^-2X(z) - C^2z^-2I(z)
>            = z^-1X(z) - Cz^-2G(z)
> where
>       G(z) = X(z)+CI(z)
> is calculated in the first line (temp0).

Correct the follow-up sign error:

(1)  I(z) = z^-1X(z) - Cz^-2G(z),
(2)  G(z) = X(z) - CI(z).

Solve (2) for X(z), substitute in (1) and simplify to get

  zI(z) = G(z),

tantamount to the third line "d = t;" in my equivalent loop. The step 
from my code to yours is called unrolling and reduces looping 
overhead. It is unrelated to the decimation that filtering only an 
input subsequence effects.


Martin

-- 
Quidquid latine dictum sit, altum viditur.

Dear all,

I am reading a piece of C code supposedly implementing a 1st order all
pass filter of transfer function A(z) = (C+z^-1)/(1+C*z^-1) with the
output decimated by a factor of 2, where C* denotes the conjugate of C*.
From the code I guess C is real, so we ignore the conjugate at the moment.
The piece of code is as follows:

     temp0 = sub(shr(in[4*i+0], 2), mult(COEFF5_1, data0));
     temp1[4*i+0] = add(data0, mult(COEFF5_1, temp0));
     data0 = sub(shr(in[4*i+2], 2), mult(COEFF5_1, temp0));
     temp1[4*i+2] = add(temp0, mult(COEFF5_1, data0));

where sub is subtraction, add is addition, mult is multiplication, shr is
right shift (can be ignored), COEFF5_1 is the coefficient C. temp1
supposedly stored the output sample of the all pass filter at 2 different
time instances. data0 stores the "filter memory" which I can't figure out
what it is.

I made the following derivation:
     Y(z) = A(z)X(z).
Since we want to do decimation, I multiply the numerator and denumerator
by (1-Cz^-1). Note that the conjugate on C is ignored. We get
     Y(z) = (1 - Cz^-1)(C + z^-1) / (1 - C^2z^-2)X(z).
Moving the denumerator to the left and then moving z^-2Y(z) term to the
right, it follows
     Y(z) = CX(z) + (1 - C^2)z^-1X(z) - Cz^-2X(z) - C^2z^-2Y(z).
To avoid buffering Y(z) (as you can see, it's named temp1 in the code,
it's just an intermidiate result), I group all terms with delay together,
call it Q(z)
     Q(z) = (1 - C^2)z^-1X(z) - Cz^-2X(z) - C^2z^-2Y(z).
Let
     I(z) = Q(z) / (1 - C^2),
we have
     Y(z) = I(z)+C[X(z)-CI(z)],
which is the first 2 lines of the code, if data0=I(z).

But then the 3rd and 4th lines are very odd. For updating I(z) from
z^-2I(z), substitute Y(z) as in the above equation into the definition of
I(z), one can get
      I(z) = z^-1X(z) - Cz^-2X(z) - C^2z^-2I(z)
           = z^-1X(z) - Cz^-2G(z)
where
      G(z) = X(z)+CI(z)
is calculated in the first line (temp0).

Does any1 have an idea how that piece of code up there implement A(z) with
decimation of factor 2? Or point to some good ref on this for me?

Many thanks.

Jiaquan


		
This message was sent using the Comp.DSP web interface on
www.DSPRelated.com