nx_av1_txfm.nx source

1// nx_av1_txfm.nx -- AV1/AV2 transform primitives. 2// 3// Layer six, the last one before actual residual reconstruction. Everything 4// here is integer and bit-exact by mandate: an inverse transform that is off 5// by one anywhere drifts the reference frames apart from the encoder's, and 6// the picture degrades over a GOP rather than failing. 7// 8// Round2 AND Round2Signed BOTH EXIST AND DIFFER. Round2 adds half and 9// arithmetic-shifts, so it rounds negatives toward positive infinity: 10// Round2(-3, 1) is -1. Round2Signed rounds AWAY from zero by construction: 11// Round2Signed(-3, 1) is -2. The spec uses each in specific places and they 12// are not interchangeable -- substituting one for the other produces a 13// transform that is correct on non-negative input and half-wrong on real 14// signal, which is exactly the kind of defect that survives a smoke test. 15// 16// THE COSINE TABLE IS 12-BIT AND QUARTER-PERIOD. Cos128_Lookup holds only 17// angles 0..64; the other three quadrants are reflections with sign flips. 18// Getting a quadrant sign wrong inverts half the basis functions, which 19// decodes to a recognisable but wrong image -- it does not look like noise, 20// so it is easy to ship. 21// 22// THE WHT IS EXACTLY INVERTIBLE. The lossless path's Walsh-Hadamard pair is 23// integer-exact in both directions, which makes it the one transform provable 24// by pure round-trip rather than against reference vectors. 25// 26// genealogy_id: av1_spec_7_13_inverse_transform 27// lineage_id: nx_av1_txfm_v1 28// license_tier: ORIGINAL 29 30import "nx_syscalls.nx" 31 32const NX_TX_COS_BITS: i64 = 12 33const NX_TX_COS_ONE: i64 = 4096 34 35// ===== rounding =================================================== 36 37func nx_av1_round2(x: i64, n: i64) -> i64 { 38 if n <= 0 { return x } 39 return (x + (1 << (n - 1))) >> n 40} 41 42func nx_av1_round2_signed(x: i64, n: i64) -> i64 { 43 if x >= 0 { return nx_av1_round2(x, n) } 44 return 0 - nx_av1_round2(0 - x, n) 45} 46 47// ===== bit reversal =============================================== 48// 49// Used to index the DCT butterfly stages. 50 51func nx_av1_brev(num_bits: i64, x: i64) -> i64 { 52 if num_bits <= 0 { return 0 } 53 var t: i64 = 0 54 var i: i64 = 0 55 while i < num_bits { 56 let bit: i64 = (x >> i) & 1 57 t = t | (bit << (num_bits - 1 - i)) 58 i = i + 1 59 } 60 return t 61} 62 63// ===== clamping =================================================== 64 65func nx_av1_clamp(v: i64, lo: i64, hi: i64) -> i64 { 66 if v < lo { return lo } 67 if v > hi { return hi } 68 return v 69} 70 71// ===== the quarter-period cosine table ============================ 72// 73// 65 entries at 12-bit precision: Cos128_Lookup[0] is 4096 (one) and 74// Cos128_Lookup[64] is 0. 75 76func nx_av1_cos_table(t: *i64) -> i64 { 77 t[0]=4096; t[1]=4095; t[2]=4091; t[3]=4085 78 t[4]=4076; t[5]=4065; t[6]=4052; t[7]=4036 79 t[8]=4017; t[9]=3996; t[10]=3973; t[11]=3948 80 t[12]=3920; t[13]=3889; t[14]=3857; t[15]=3822 81 t[16]=3784; t[17]=3745; t[18]=3703; t[19]=3659 82 t[20]=3612; t[21]=3564; t[22]=3513; t[23]=3461 83 t[24]=3406; t[25]=3349; t[26]=3290; t[27]=3229 84 t[28]=3166; t[29]=3102; t[30]=3035; t[31]=2967 85 t[32]=2896; t[33]=2824; t[34]=2751; t[35]=2675 86 t[36]=2598; t[37]=2520; t[38]=2440; t[39]=2359 87 t[40]=2276; t[41]=2191; t[42]=2106; t[43]=2019 88 t[44]=1931; t[45]=1842; t[46]=1751; t[47]=1660 89 t[48]=1567; t[49]=1474; t[50]=1380; t[51]=1285 90 t[52]=1189; t[53]=1092; t[54]=995; t[55]=897 91 t[56]=799; t[57]=700; t[58]=601; t[59]=501 92 t[60]=401; t[61]=301; t[62]=201; t[63]=101 93 t[64]=0 94 return 65 95} 96 97// cos128 over the full 256-step circle, by quadrant reflection 98func nx_av1_cos128(t: *i64, angle: i64) -> i64 { 99 let a: i64 = angle & 255 100 if a <= 64 { return t[a] } 101 if a <= 128 { return 0 - t[128 - a] } 102 if a <= 192 { return 0 - t[a - 128] } 103 return t[256 - a] 104} 105 106func nx_av1_sin128(t: *i64, angle: i64) -> i64 { 107 return nx_av1_cos128(t, angle - 64) 108} 109 110// ===== the lossless Walsh-Hadamard pair =========================== 111// 112// Exactly invertible in integer arithmetic. The forward scales by 4 and the 113// inverse shifts it back, so the pair composes to identity. 114 115func nx_av1_fwht4(inp: *i64, out: *i64) -> i64 { 116 var a: i64 = inp[0] 117 var b: i64 = inp[1] 118 var c: i64 = inp[2] 119 var d: i64 = inp[3] 120 a = a + b 121 d = d - c 122 let e: i64 = (a - d) >> 1 123 b = e - b 124 c = e - c 125 a = a - c 126 d = d + b 127 out[0] = a 128 out[1] = c 129 out[2] = d 130 out[3] = b 131 return 1 132} 133 134func nx_av1_iwht4(inp: *i64, out: *i64, shift: i64) -> i64 { 135 var a: i64 = inp[0] >> shift 136 var c: i64 = inp[1] >> shift 137 var d: i64 = inp[2] >> shift 138 var b: i64 = inp[3] >> shift 139 a = a + c 140 d = d - b 141 let e: i64 = (a - d) >> 1 142 b = e - b 143 c = e - c 144 a = a - b 145 d = d + c 146 out[0] = a 147 out[1] = b 148 out[2] = c 149 out[3] = d 150 return 1 151} 152 153// ===== identity transform scale factors =========================== 154// 155// Each size has its OWN factor and they are not powers of two throughout: 156// 4 and 16 use the 5793/4096 root-two scaling, 8 and 32 are plain doublings. 157// Applying the wrong one leaves the residual off by a constant ratio. 158 159func nx_av1_identity4(v: i64) -> i64 { 160 return nx_av1_round2(v * 5793, NX_TX_COS_BITS) 161} 162 163func nx_av1_identity8(v: i64) -> i64 { 164 return v * 2 165} 166 167func nx_av1_identity16(v: i64) -> i64 { 168 return nx_av1_round2(v * 2 * 5793, NX_TX_COS_BITS) 169} 170 171func nx_av1_identity32(v: i64) -> i64 { 172 return v * 4 173} 174 175// ===== the butterfly primitives =================================== 176// 177// B() rotates a pair by an angle; H() is the plain Hadamard add/subtract. 178// Both write back through the same array, which is how the spec expresses 179// the in-place butterfly network. 180 181func nx_av1_butterfly_b(t: *i64, arr: *i64, a: i64, b: i64, angle: i64, flip: i64) -> i64 { 182 let x: i64 = arr[a] 183 let y: i64 = arr[b] 184 let c: i64 = nx_av1_cos128(t, angle) 185 let s: i64 = nx_av1_sin128(t, angle) 186 if flip == 1 { 187 arr[a] = nx_av1_round2(x * s + y * c, NX_TX_COS_BITS) 188 arr[b] = nx_av1_round2(x * c - y * s, NX_TX_COS_BITS) 189 } else { 190 arr[a] = nx_av1_round2(x * c - y * s, NX_TX_COS_BITS) 191 arr[b] = nx_av1_round2(x * s + y * c, NX_TX_COS_BITS) 192 } 193 return 1 194} 195 196func nx_av1_butterfly_h(arr: *i64, a: i64, b: i64, flip: i64) -> i64 { 197 var i: i64 = a 198 var j: i64 = b 199 if flip == 1 { i = b; j = a } 200 let x: i64 = arr[i] 201 let y: i64 = arr[j] 202 arr[i] = x + y 203 arr[j] = x - y 204 return 1 205}