Search tips

# Discussion Groups | Comp.DSP | 1st order all pass filter implementation

There are 2 messages in this thread.

You are currently looking at messages 1 to .

Is this discussion worth a thumbs up?

0

# 1st order all pass filter implementation - Huo Jiaquan - 2005-07-19 07:56:00

```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));
data0 = sub(shr(in[4*i+2], 2), mult(COEFF5_1, temp0));

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.

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```

# Re: 1st order all pass filter implementation - Martin Eisenberg - 2005-07-24 12:11:00

```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.```