diff options
Diffstat (limited to 'appl/wm/mpeg/fixidct.b')
| -rw-r--r-- | appl/wm/mpeg/fixidct.b | 188 |
1 files changed, 188 insertions, 0 deletions
diff --git a/appl/wm/mpeg/fixidct.b b/appl/wm/mpeg/fixidct.b new file mode 100644 index 00000000..992cf837 --- /dev/null +++ b/appl/wm/mpeg/fixidct.b @@ -0,0 +1,188 @@ +implement IDCT; + +include "sys.m"; +include "mpegio.m"; + +init() +{ +} + +# IDCT based on Arai, Agui, and Nakajima, using flow chart Figure 4.8 +# of Pennebaker & Mitchell, JPEG: Still Image Data Compression Standard. +# Remember IDCT is reverse of flow of DCT. +# Nasty truncated integer version (not compliant). + +B0: con 16; +B1: con 16; +M: con (1 << B0); +N: con (1 << B1); + +a0: con 1.414; +a1: con 0.707; +a2: con 0.541; +a3: con 0.707; +a4: con 1.307; +a5: con -0.383; + +A0: con int (a0 * real N); +A1: con int (a1 * real M); +A2: con int (a2 * real M); +A3: con int (a3 * real M); +A4: con int (a4 * real M); +A5: con int (a5 * real M); + +# scaling factors from eqn 4-35 of P&M +s1: con 1.0196; +s2: con 1.0823; +s3: con 1.2026; +s4: con 1.4142; +s5: con 1.8000; +s6: con 2.6131; +s7: con 5.1258; + +S1: con int (s1 * real N); +S2: con int (s2 * real N); +S3: con int (s3 * real N); +S4: con int (s4 * real N); +S5: con int (s5 * real N); +S6: con int (s6 * real N); +S7: con int (s7 * real N); + +# overall normalization of 1/16, folded into premultiplication on vertical pass +S: con 4; +scale: con 0.0625; + +idct(b: array of int) +{ + x, y: int; + + r := array[8*8] of int; + + # transform horizontally + for(y=0; y<8; y++){ + eighty := y<<3; + # if all non-DC components are zero, just propagate the DC term + if(b[eighty+1]==0) + if(b[eighty+2]==0 && b[eighty+3]==0) + if(b[eighty+4]==0 && b[eighty+5]==0) + if(b[eighty+6]==0 && b[eighty+7]==0){ + v := b[eighty]*A0; + r[eighty+0] = v; + r[eighty+1] = v; + r[eighty+2] = v; + r[eighty+3] = v; + r[eighty+4] = v; + r[eighty+5] = v; + r[eighty+6] = v; + r[eighty+7] = v; + continue; + } + + # step 5 + in1 := S1*b[eighty+1]; + in3 := S3*b[eighty+3]; + in5 := S5*b[eighty+5]; + in7 := S7*b[eighty+7]; + f2 := S2*b[eighty+2]; + f3 := S6*b[eighty+6]; + f5 := (in1+in7); + f7 := (in5+in3); + + # step 4 + g2 := f2-f3; + g4 := (in5-in3); + g6 := (in1-in7); + g7 := f5+f7; + + # step 3.5 + t := ((g4+g6)>>B0)*A5; + + # step 3 + f0 := A0*b[eighty+0]; + f1 := S4*b[eighty+4]; + f3 += f2; + f2 = A1*(g2>>B0); + + # step 2 + g0 := f0+f1; + g1 := f0-f1; + g3 := f2+f3; + g4 = t-A2*(g4>>B0); + g5 := A3*((f5-f7)>>B0); + g6 = A4*(g6>>B0)+t; + + # step 1 + f0 = g0+g3; + f1 = g1+f2; + f2 = g1-f2; + f3 = g0-g3; + f5 = g5-g4; + f6 := g5+g6; + f7 = g6+g7; + + # step 6 + r[eighty+0] = (f0+f7); + r[eighty+1] = (f1+f6); + r[eighty+2] = (f2+f5); + r[eighty+3] = (f3-g4); + r[eighty+4] = (f3+g4); + r[eighty+5] = (f2-f5); + r[eighty+6] = (f1-f6); + r[eighty+7] = (f0-f7); + } + + # transform vertically + for(x=0; x<8; x++){ + # step 5 + in1 := S1*(r[x+8]>>(B1+S)); + in3 := S3*(r[x+24]>>(B1+S)); + in5 := S5*(r[x+40]>>(B1+S)); + in7 := S7*(r[x+56]>>(B1+S)); + f2 := S2*(r[x+16]>>(B1+S)); + f3 := S6*(r[x+48]>>(B1+S)); + f5 := (in1+in7); + f7 := (in5+in3); + + # step 4 + g2 := f2-f3; + g4 := (in5-in3); + g6 := (in1-in7); + g7 := f5+f7; + + # step 3.5 + t := ((g4+g6)>>B0)*A5; + + # step 3 + f0 := A0*(r[x]>>(B1+S)); + f1 := S4*(r[x+32]>>(B1+S)); + f3 += f2; + f2 = A1*(g2>>B0); + + # step 2 + g0 := f0+f1; + g1 := f0-f1; + g3 := f2+f3; + g4 = t-A2*(g4>>B0); + g5 := A3*((f5-f7)>>B0); + g6 = A4*(g6>>B0)+t; + + # step 1 + f0 = g0+g3; + f1 = g1+f2; + f2 = g1-f2; + f3 = g0-g3; + f5 = g5-g4; + f6 := g5+g6; + f7 = g6+g7; + + # step 6 + b[x] = (f0+f7)>>B1; + b[x+8] = (f1+f6)>>B1; + b[x+16] = (f2+f5)>>B1; + b[x+24] = (f3-g4)>>B1; + b[x+32] = (f3+g4)>>B1; + b[x+40] = (f2-f5)>>B1; + b[x+48] = (f1-f6)>>B1; + b[x+56] = (f0-f7)>>B1; + } +} |