Martin Eisenberg
07-24-2005, 05:11 PM
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.
> 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.