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

# 1st order all pass filter implementation

Started by ●July 19, 2005

Posted by ●July 24, 2005

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 2More 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.

Posted by ●July 24, 2005

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 2More 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.

Posted by ●July 19, 2005

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