1# mach: bfin 2 3// GENERIC BIQUAD: 4// --------------- 5// x ---------+---------|---------+-------y 6// | |t1 | 7// | D | 8// | a1 | b1 | 9// +---<-----|---->----+ 10// | | | 11// | D | D's are delays 12// | a2 | b2 | ">" represent multiplications 13// +---<-----|---->----+ 14// To test this routine, use a biquad with a pole pair at z = (0.7 +- 0.1j), 15// and a double zero at z = -1.0, which is a low-pass. The transfer function is: 16// 1 + 2z^-1 + z^-2 17// H(z) = ---------------------- 18// 1 - 1.4z^-1 + 0.5z^-2 19// a1 = 1.4 20// a2 = -0.5 21// b1 = 2 22// b2 = 1 23// This filter conforms to the biquad test in BDT, since it has coefficients 24// larger than 1.0 in magnitude, and b0=1. (Note that the a's have a negative 25// sign.) 26// This filter can be simulated in matlab. To simulate one biquad, use 27// A = [1.0, -1.4, 0.5] 28// B = [1, 2, 1] 29// Y=filter(B,A,X) 30// To simulate two cascaded biquads, use 31// Y=filter(B,A,filter(B,A,X)) 32// SCALED COEFFICIENTS: 33// -------------------- 34// In order to conform to 1.15 representation, must scale coeffs by 0.5. 35// This requires an additional internal re-scale. The equations for the Type II 36// filter are: 37// t1 = x + a1*t1*z^-1 + a2*t1*z^-2 38// y = b0*t1 + b1*t1*z^-1 + b2*t1*z^-2 39// (Note inclusion of term b0, which in the example is b0 = 1.) 40// If all coeffs are replaced by 41// ai --> ai' = 0.5*a1 42// then the two equations become 43// t1 = x + 2*a1'*t1*z^-1 + 2*a2'*t1*z^-2 44// 0.5*y = b0'*t1 + b1'*t1*z^-1 + b2'*t1*z^-2 45// which can be implemented as: 46// 2.0 b0'=0.5 47// x ---------+--->-----|---->----+-------y 48// | |t1 | 49// | D | 50// | a1' | b1' | 51// +---<-----|---->----+ 52// | | | 53// | D | 54// | a2' | b2' | 55// +---<-----|---->----+ 56// But, b0' can be eliminated by: 57// x ---------+---------|---------+-------y 58// | | | 59// | V 2.0 | 60// | | | 61// | |t1 | 62// | D | 63// | a1' | b1' | 64// +---<-----|---->----+ 65// | | | 66// | D | 67// | a2' | b2' | 68// +---<-----|---->----+ 69// Function biquadf() computes this implementation on float data. 70// CASCADED BIQUADS 71// ---------------- 72// Cascaded biquads are simulated by simply cascading copies of the 73// filter defined above. However, one must be careful with the resulting 74// filter, as it is not very stable numerically (double poles in the 75// vecinity of +1). It would of course be better to cascade different 76// filters, as that would result in more stable structures. 77// The functions biquadf() and biquadR() have been tested with up to 3 78// stages using this technique, with inputs having small signal amplitude 79// (less than 0.001) and under 300 samples. 80// 81// In order to pipeline, need to maintain two pointers into the state 82// array: one to load (I0) and one to store (I2). This is required since 83// the load of iteration i+1 is hoisted above the store of iteration i. 84 85.include "testutils.inc" 86 start 87 88 89 // I3 points to input buffer 90 loadsym I3, input; 91 92 // P1 points to output buffer 93 loadsym P1, output; 94 95 R0 = 0; R7 = 0; 96 97 P2 = 10; 98 LSETUP ( L$0 , L$0end ) LC0 = P2; 99L$0: 100 101 // I0 and I2 are pointers to state 102 loadsym I0, state; 103 I2 = I0; 104 105 // pointer to coeffs 106 loadsym I1, Coeff; 107 108 R0.H = W [ I3 ++ ]; // load input value into RH0 109 A0.w = R0; // A0 holds x 110 111 P2 = 2; 112 LSETUP ( L$1 , L$1end ) LC1 = P2; 113 114 // load 2 coeffs into R1 and R2 115 // load state into R3 116 R1 = [ I1 ++ ]; 117 MNOP || R2 = [ I1 ++ ] || R3 = [ I0 ++ ]; 118 119L$1: 120 121 // A1=b1*s0 A0=a1*s0+x 122 A1 = R1.L * R3.L, A0 += R1.H * R3.L || R1 = [ I1 ++ ] || NOP; 123 124 // A1+=b2*s1 A0+=a2*s1 125 // and move scaled value in A0 (t1) into RL4 126 A1 += R2.L * R3.H, R4.L = ( A0 += R2.H * R3.H ) (S2RND) || R2 = [ I1 ++ ] || NOP; 127 128 // Advance state. before: 129 // R4 = uuuu t1 130 // R3 = stat[1] stat[0] 131 // after PACKLL: 132 // R3 = stat[0] t1 133 R5 = PACK( R3.L , R4.L ) || R3 = [ I0 ++ ] || NOP; 134 135 // collect output into A0, and move to RL0. 136 // Keep output value in A0, since it is also 137 // the accumulator used to store the input to 138 // the next stage. Also, store updated state 139L$1end: 140 R0.L = ( A0 += A1 ) || [ I2 ++ ] = R5 || NOP; 141 142 // store output 143L$0end: 144 W [ P1 ++ ] = R0; 145 146 // Check results 147 loadsym I2, output; 148 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x0028 ); 149 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x0110 ); 150 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x0373 ); 151 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x075b ); 152 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x0c00 ); 153 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x1064 ); 154 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x13d3 ); 155 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x15f2 ); 156 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x16b9 ); 157 R0.L = W [ I2 ++ ]; DBGA ( R0.L , 0x1650 ); 158 159 pass 160 161 .data 162state: 163 .dw 0x0000 164 .dw 0x0000 165 .dw 0x0000 166 .dw 0x0000 167 .dw 0x0000 168 .dw 0x0000 169 .dw 0x0000 170 .dw 0x0000 171 172 .data 173Coeff: 174 .dw 0x7fff 175 .dw 0x5999 176 .dw 0x4000 177 .dw 0xe000 178 .dw 0x7fff 179 .dw 0x5999 180 .dw 0x4000 181 .dw 0xe000 182input: 183 .dw 0x0028 184 .dw 0x0000 185 .dw 0x0000 186 .dw 0x0000 187 .dw 0x0000 188 .dw 0x0000 189 .dw 0x0000 190 .dw 0x0000 191 .dw 0x0000 192 .dw 0x0000 193 .dw 0x0000 194output: 195 .dw 0x0000 196 .dw 0x0000 197 .dw 0x0000 198 .dw 0x0000 199 .dw 0x0000 200 .dw 0x0000 201 .dw 0x0000 202 .dw 0x0000 203 .dw 0x0000 204 .dw 0x0000 205 .dw 0x0000 206 .dw 0x0000 207 .dw 0x0000 208