8e2be5d59fa2fcbe34d0ca310424f5729643683e
[daala.git] / doc / ietf / draft-terriberry-netvc-codingtools.xml
1 <?xml version="1.0" encoding="utf-8"?>
2 <!DOCTYPE rfc SYSTEM 'rfc2629.dtd'>
3 <?rfc toc="yes" symrefs="yes" ?>
4
5 <rfc ipr="trust200902" category="info" docName="draft-terriberry-netvc-codingtools-01">
6
7 <front>
8 <title abbrev="Coding Tools">Coding Tools for a Next Generation Video Codec</title>
9 <author initials="T.B." surname="Terriberry" fullname="Timothy B. Terriberry">
10 <organization>Mozilla Corporation</organization>
11 <address>
12 <postal>
13 <street>331 E. Evelyn Avenue</street>
14 <city>Mountain View</city>
15 <region>CA</region>
16 <code>94041</code>
17 <country>USA</country>
18 </postal>
19 <phone>+1 650 903-0800</phone>
20 <email>tterribe@xiph.org</email>
21 </address>
22 </author>
23 <author initials="N.E." surname="Egge" fullname="Nathan E. Egge">
24 <organization>Mozilla Corporation</organization>
25 <address>
26 <postal>
27 <street>331 E. Evelyn Avenue</street>
28 <city>Mountain View</city>
29 <region>CA</region>
30 <code>94041</code>
31 <country>USA</country>
32 </postal>
33 <phone>+1 650 903-0800</phone>
34 <email>negge@xiph.org</email>
35 </address>
36 </author>
37
38 <date day="31" month="October" year="2016"/>
39 <area>ART</area>
40 <workgroup>netvc</workgroup>
41
42 <abstract>
43 <t>
44 This document proposes a number of coding tools that could be incorporated into
45  a next-generation video codec.
46 </t>
47 </abstract>
48 </front>
49
50 <middle>
51 <section anchor="intro" title="Introduction">
52 <t>
53 One of the biggest contributing factors to the success of the Internet is that
54  the underlying protocols are implementable on a royalty-free basis.
55 This allows them to be implemented widely and easily distributed by application
56  developers, service operators, and end users, without asking for permission.
57 In order to produce a next-generation video codec that is competitive with the
58  best patent-encumbered standards, yet avoids patents which are not available
59  on an open-source compatible, royalty-free basis, we must use old coding tools
60  in new ways and develop new coding tools.
61 This draft documents some of the tools we have been working on for inclusion in
62  such a codec.
63 This is early work, and the performance of some of these tools (especially in
64  relation to other approaches) is not yet fully known.
65 Nevertheless, it still serves to outline some possibilities an eventual working
66  group, if formed, could consider.
67 </t>
68
69 </section>
70
71 <section anchor="entropy_coding" title="Entropy Coding">
72 <t>
73 The basic theory of entropy coding was well-established by the late
74  1970's&nbsp;<xref target="Pas76"/>.
75 Modern video codecs have focused on Huffman codes (or "Variable-Length
76  Codes"/VLCs) and binary arithmetic coding.
77 Huffman codes are limited in the amount of compression they can provide and the
78  design flexibility they allow, but as each code word consists of an integer
79  number of bits, their implementation complexity is very low, so they were
80  provided at least as an option in every video codec up through H.264.
81 Arithmetic coding, on the other hand, uses code words that can take up
82  fractional parts of a bit, and are more complex to implement.
83 However, the prevalence of cheap, H.264 High Profile hardware, which requires
84  support for arithmetic coding, shows that it is no longer so expensive that a
85  fallback VLC-based approach is required.
86 Having a single entropy-coding method simplifies both up-front design costs and
87  interoperability.
88 </t>
89 <t>
90 However, the primary limitation of arithmetic coding is that it is an
91  inherently serial operation.
92 A given symbol cannot be decoded until the previous symbol is decoded, because
93  the bits (if any) that are output depend on the exact state of the decoder at
94  the time it is decoded.
95 This means that a hardware implementation must run at a sufficiently high clock
96  rate to be able to decode all of the symbols in a frame.
97 Higher clock rates lead to increased power consumption, and in some cases the
98  entropy coding is actually becoming the limiting factor in these designs.
99 </t>
100 <t>
101 As fabrication processes improve, implementers are very willing to trade
102  increased gate count for lower clock speeds.
103 So far, most approaches to allowing parallel entropy coding have focused on
104  splitting the encoded symbols into multiple streams that can be decoded
105  independently.
106 This "independence" requirement has a non-negligible impact on compression,
107  parallelizability, or both.
108 For example, H.264 can split frames into "slices" which might cover only a
109  small subset of the blocks in the frame.
110 In order to allow decoding these slices independently, they cannot use context
111  information from blocks in other slices (harming compression).
112 Those contexts must adapt rapidly to account for the generally small number of
113  symbols available for learning probabilities (also harming compression).
114 In some cases the number of contexts must be reduced to ensure enough symbols
115  are coded in each context to usefully learn probabilities at all (once more,
116  harming compression).
117 Furthermore, an encoder must specially format the stream to use multiple slices
118  per frame to allow any parallel entropy decoding at all.
119 Encoders rarely have enough information to evaluate this "compression
120  efficiency" vs. "parallelizability" trade-off, since they don't generally know
121  the limitations of the decoders for which they are encoding.
122 That means there will be many files or streams which could have been decoded if
123  they were encoded with different options, but which a given decoder cannot
124  decode because of bad choices made by the encoder (at least from the
125  perspective of that decoder).
126 The same set of drawbacks apply to the DCT token partitions in
127  VP8&nbsp;<xref target="RFC6386"/>.
128 </t>
129
130 <section anchor="nonbinary_coding" title="Non-binary Arithmetic Coding">
131 <t>
132 Instead, we propose a very different approach: use non-binary arithmetic
133  coding.
134 In binary arithmetic coding, each decoded symbol has one of two possible
135  values: 0 or 1.
136 The original arithmetic coding algorithms allow a symbol to take on any number
137  of possible values, and allow the size of that alphabet to change with each
138  symbol coded.
139 Reasonable values of N (for example, N&nbsp;&lt;=&nbsp;16) offer the potential
140  for a decent throughput increase for a reasonable increase in gate count for
141  hardware implementations.
142 </t>
143 <t>
144 Binary coding allows a number of computational simplifications.
145 For example, for each coded symbol, the set of valid code points is partitioned
146  in two, and the decoded value is determined by finding the partition in which
147  the actual code point that was received lies.
148 This can be determined by computing a single partition value (in both the
149  encoder and decoder) and (in the decoder) doing a single comparison.
150 A non-binary arithmetic coder partitions the set of valid code points
151  into multiple pieces (one for each possible value of the coded symbol).
152 This requires the encoder to compute two partition values, in general (for both
153  the upper and lower bound of the symbol to encode).
154 The decoder, on the other hand, must search the partitions for the one that
155  contains the received code point.
156 This requires computing at least O(log&nbsp;N) partition values.
157 </t>
158 <t>
159 However, coding a parameter with N possible values with a binary arithmetic
160  coder requires O(log&nbsp;N) symbols in the worst case (the only case
161  that matters for hardware design).
162 Hence, this does not represent any actual savings (indeed, it represents an
163  increase in the number of partition values computed by the encoder).
164 In addition, there are a number of overheads that are per-symbol, rather than
165  per-value.
166 For example, renormalization (which enlarges the set of valid code points after
167  partitioning has reduced it too much), carry propagation (to deal with the
168  case where the high and low ends of a partition straddle a bit boundary),
169  etc., are all performed on a symbol-by-symbol basis.
170 Since a non-binary arithmetic coder codes a given set of values with fewer
171  symbols than a binary one, it incurs these per-symbol overheads less often.
172 This suggests that a non-binary arithmetic coder can actually be more efficient
173  than a binary one.
174 </t>
175
176 </section>
177
178 <section anchor="nonbinary_modeling" title="Non-binary Context Modeling">
179 <t>
180 The other aspect that binary coding simplifies is probability modeling.
181 In arithmetic coding, the size of the sets the code points are partitioned into
182  are (roughly) proportional to the probability of each possible symbol value.
183 Estimating these probabilities is part of the coding process, though it can be
184  cleanly separated from the task of actually producing the coded bits.
185 In a binary arithmetic coder, this requires estimating the probability of only
186  one of the two possible values (since the total probability is 1.0).
187 This is often done with a simple table lookup that maps the old probability and
188  the most recently decoded symbol to a new probability to use for the next
189  symbol in the current context.
190 The trade-off, of course, is that non-binary symbols must be "binarized" into
191  a series of bits, and a context (with an associated probability) chosen for
192  each one.
193 </t>
194 <t>
195 In a non-binary arithmetic coder, the decoder must compute at least
196  O(log&nbsp;N) cumulative probabilities (one for each partition value it
197  needs).
198 Because these probabilities are usually not estimated directly in "cumulative"
199  form, this can require computing (N&nbsp;-&nbsp;1) non-cumulative probability
200  values.
201 Unless N is very small, these cannot be updated with a single table lookup.
202 The normal approach is to use "frequency counts".
203 Define the frequency of value k to be
204 <figure align="center">
205 <artwork align="center"><![CDATA[
206 f[k] = A*<the number of times k has been observed> + B
207 ]]></artwork>
208 </figure>
209  where A and B are parameters (usually A=2 and B=1 for a traditional
210  Krichevsky-Trofimov estimator).
211 The resulting probability, p[k], is given by
212 <figure align="center">
213 <artwork align="center"><![CDATA[
214        N-1
215        __
216   ft = \   f[k]
217        /_
218        k=0
219
220        f[k]
221 p[k] = ----
222         ft
223 ]]></artwork>
224 </figure>
225 When ft grows too large, the frequencies are rescaled (e.g., halved, rounding
226  up to prevent reduction of a probability to 0).
227 </t>
228 <t>
229 When ft is not a power of two, partitioning the code points requires actual
230  divisions (see <xref target="RFC6716"/> Section&nbsp;4.1 for one detailed
231  example of exactly how this is done).
232 These divisions are acceptable in an audio codec like
233  Opus&nbsp;<xref target="RFC6716"/>, which only has to code a few hundreds of
234  these symbols per second.
235 But video requires hundreds of thousands of symbols per second, at a minimum,
236  and divisions are still very expensive to implement in hardware.
237 </t>
238 <t>
239 There are two possible approaches to this.
240 One is to come up with a replacement for frequency counts that produces
241  probabilities that sum to a power of two.
242 Some possibilities, which can be applied individually or in combination:
243 <list style="numbers">
244 <t>
245 Use probabilities that are fixed for the duration of a frame.
246 This is the approach taken by VP8, for example, even though it uses a binary
247  arithmetic coder.
248 In fact, it is possible to convert many of VP8's existing binary-alphabet
249  probabilities into probabilities for non-binary alphabets, an approach that is
250  used in the experiment presented at the end of this section.
251 </t>
252 <t>
253 Use parametric distributions.
254 For example, DCT coefficient magnitudes usually have an approximately
255  exponential distribution.
256 This distribution can be characterized by a single parameter, e.g., the
257  expected value.
258 The expected value is trivial to update after decoding a coefficient.
259 For example
260 <figure align="center">
261 <artwork align="center"><![CDATA[
262 E[x[n+1]] = E[x[n]] + floor(C*(x[n] - E[x[n]]))
263 ]]></artwork>
264 </figure>
265  produces an exponential moving average with a decay factor of
266  (1&nbsp;-&nbsp;C).
267 For a choice of C that is a negative power of two (e.g., 1/16 or 1/32 or
268  similar), this can be implemented with two adds and a shift.
269 Given this expected value, the actual distribution to use can be obtained from
270  a small set of pre-computed distributions via a lookup table.
271 Linear interpolation between these pre-computed values can improve accuracy, at
272  the cost of O(N) computations, but if N is kept small this is trivially
273  parallelizable, in SIMD or otherwise.
274 </t>
275 <t>
276 Change the frequency count update mechanism so that ft is constant.
277 This approach is described in the next section.
278 </t>
279 </list>
280 </t>
281 </section>
282 <section anchor="dyadic_adaptation" title="Dyadic Adaptation">
283 <t>
284 The goal with context adaptation using dyadic probabilities is to maintain
285  the invariant that the probabilities all sum to a power of two before and
286  after adaptation.
287 This can be achieved with a special update function that blends the cumulative
288  probabilities of the current context with a cumulative distribution function
289  where the coded symbol has probability 1.
290 </t>
291 <t>
292 Suppose we have model for a given context that codes 8 symbols with the
293  following probabilities:
294 <figure align="center">
295 <artwork align="center"><![CDATA[
296 +------+------+------+------+------+------+------+------+
297 | p[0] | p[1] | p[2] | p[3] | p[4] | p[5] | p[6] | p[7] |
298 +------+------+------+------+------+------+------+------+
299 |  1/8 |  1/8 | 3/16 | 1/16 | 1/16 | 3/16 |  1/8 |  1/8 |
300 +------+------+------+------+------+------+------+------+
301 ]]></artwork>
302 </figure>
303 Then the cumulative distribution function is:
304 <figure align="center">
305 <artwork align="left"><![CDATA[
306    CDF
307
308  1  +                                                +------+
309     |                                                |
310     |                                         +------+
311     |                                         |
312 3/4 +                                  +------+
313     |                                  |
314     |                                  |
315     |                           +------+
316 1/2 +                    +------+
317     |             +------+
318     |             |
319     |             |
320 1/4 +      +------+
321     |      |
322     +------+
323     |
324  0  +------+------+------+------+------+------+------+------+ Bin
325      fl[1]  fl[2]  fl[3]  fl[4]  fl[5]  fl[6]  fl[7]  fl[8]
326 ]]></artwork>
327 </figure>
328 Suppose we code symbol 3 and wish to update the context model so that this
329  symbol is now more likely.
330 This can be done by blending the CDF for the current context with a CDF
331  that has symbol 3 with likelihood 1.
332 <figure align="center">
333 <artwork align="left"><![CDATA[
334    CDF
335
336  1  +                    +----------------------------------+
337     |                    |
338     |                    |
339     |                    |
340  0  +------+------+------+------+------+------+------+------+ Bin
341      fl[1]  fl[2]  fl[3]  fl[4]  fl[5]  fl[6]  fl[7]  fl[8]
342 ]]></artwork>
343 </figure>
344 Given an adaptation rate g between 0 and 1, and assuming ft = 2^4 = 16, what
345  we are computing is:
346 <figure align="center">
347 <artwork align="center"><![CDATA[
348 +------+------+------+------+------+------+------+------+
349 |   2  |   4  |   7  |   8  |   9  |  12  |  14  |  16  |  * (1 - g)
350 +------+------+------+------+------+------+------+------+
351
352                             +
353
354 +------+------+------+------+------+------+------+------+
355 |   0  |   0  |   0  |  16  |  16  |  16  |  16  |  16  |  * g
356 +------+------+------+------+------+------+------+------+
357 ]]></artwork>
358 </figure>
359 In order to prevent the probability of any one symbol from going to zero, the
360  blending functions above and below the coded symbol are adjusted so that no
361  adjacent cumulative probabilities are the same.
362 </t>
363 <t>
364 Let M be the alphabet size and 1/2^r be the adaptation rate:
365 </t>
366 <t>
367 <figure align="center">
368 <artwork align="center"><![CDATA[
369         ( fl[i] - floor((fl[i] + 2^r - i - 1)/2^r), i <= coded symbol
370 fl[i] = <
371         ( fl[i] - floor((fl[i] + M - i - ft)/2^r),  i > coded symbol
372 ]]></artwork>
373 </figure>
374 Applying these formulas to the example CDF where M = 8 with adaptation rate
375  1/2^16 gives the updated CDF:
376 <figure align="center">
377 <artwork align="center"><![CDATA[
378 +------+------+------+------+------+------+------+------+
379 |   1  |   3  |   6  |   9  |  10  |  13  |  15  |  16  |
380 +------+------+------+------+------+------+------+------+
381 ]]></artwork>
382 </figure>
383 Looking at the graph of the CDF we see that the likelihood for symbol 3
384  has gone up from 1/16 to 3/16, dropping the likelihood of all other symbols
385  to make room.
386 <figure align="center">
387 <artwork align="left"><![CDATA[
388    CDF
389
390  1  +                                                +------+
391     |                                         +------+
392     |                                         |
393     |                                  +------+
394 3/4 +                                  |
395     |                                  |
396     |                           +------+
397     |                    +------+
398 1/2 +                    |
399     |                    |
400     |             +------+
401     |             |
402 1/4 +             |
403     |      +------+
404     |      |
405     +------+
406  0  +------+------+------+------+------+------+------+------+ Bin
407      fl[1]  fl[2]  fl[3]  fl[4]  fl[5]  fl[6]  fl[7]  fl[8]
408 ]]></artwork>
409 </figure>
410 </t>
411 </section>
412
413 <section title="Simplified Partition Function">
414 <t>
415 Rather than changing the context modeling, the other approach is to change the
416  function used to partition the set of valid code points so that it does not
417  need a division, even when ft is not a power of two.
418 Let the range of valid code points in the current arithmetic coder state be
419  [L,&nbsp;L&nbsp;+&nbsp;R), where L is the lower bound of the range and R is
420  the number of valid code points.
421 Assume that ft&nbsp;&lt;=&nbsp;R&nbsp;&lt;&nbsp;2*ft (this is easy to enforce
422  with the normal rescaling operations used with frequency counts).
423 Then one possible partition function is
424 <figure align="center">
425 <artwork align="center"><![CDATA[
426 r[k] = fl[k] + min(fl[k], R - ft)
427 ]]></artwork>
428 </figure>
429  so that the new range after coding symbol k is
430  [L&nbsp;+&nbsp;r[k],&nbsp;L&nbsp;+&nbsp;r[k+1]).
431 </t>
432 <t>
433 This is a variation of the partition function proposed
434  by&nbsp;<xref target="SM98"/>.
435 The size of the new partition (r[k+1]&nbsp;-&nbsp;r[k]) is no longer truly
436  proportional to R*p[k].
437 This partition function counts values of fl[k] smaller than R&nbsp;-&nbsp;ft
438  double compared to values larger than R&nbsp;-&nbsp;ft.
439 This over-estimates the probability of symbols at the start of the alphabet
440  and underestimates the probability of symbols at the end of the alphabet.
441 The amount of the range allocated to a symbol can be off by up to a factor of
442  2 compared to its fair share, implying a peak error as large as one bit per
443  symbol.
444 However, if the probabilities are accurate and the symbols being coded are
445  independent, the average inefficiency introduced will be as low as
446  log2(log2(e)*2/e)&nbsp;~=&nbsp;0.0861 bits per symbol.
447 This error can, of course, be reduced by coding fewer symbols with larger
448  alphabets.
449 In practice the overhead is roughly equal to the overhead introduced by other
450  approximate arithmetic coders like H.264's CABAC.
451 However, probabilities near one-half tend to have the most overhead.
452 In fact, probabilities in the range of 40% to 60% for a binary symbol may not
453  be worth modeling, since the compression gains may be entirely countered
454  by the added overhead, making it cheaper and faster to code such values as
455  raw bits.
456 This problem is partially alleviated by using larger alphabets.
457 </t>
458 <section title="Reduced-Overhead Partition Function">
459 <t>
460 A slightly more complicated partition function can reduce the overhead while
461  still avoiding the division.
462 This is done by splitting things into two cases:
463 <list>
464 <t>
465 Case 1: <![CDATA[R - ft > ft - (R - ft)]]><vspace blankLines="1"/>
466 That is, we have more values that are double-counted than single-counted.
467 In this case, we still double-count the first 2*R&nbsp;-&nbsp;3*ft values,
468  but after that we alternate between single-counting and double-counting
469  for the rest.
470 </t>
471 <t>
472 Case 2: <![CDATA[R - ft < ft - (R - ft)]]><vspace blankLines="1"/>
473 That is, we have more values that are single-counted than double-counted.
474 In this case, we alternate between single-counting and double-counting for
475  the first 2*(R&nbsp;-&nbsp;ft) values, and single-count the rest.
476 </t>
477 </list>
478 For two equiprobable symbols in different places in the alphabet, this
479  reduces the maximum ratio of over-estimation to under-estimation from 2:1
480  for the previous partition function to either 4:3 or 3:2 (for each of the
481  two cases above, respectively), assuming symbol probabilities significantly
482  greater than the minimum possible.
483 That reduces the worst-case per-symbol overhead from 1 bit to 0.58 bits.
484 </t>
485 <t>
486 The resulting reduced-overhead partition function is
487 </t>
488 <figure align="center">
489 <artwork align="center"><![CDATA[
490    e = max(2*R - 3*ft, 0)
491 r[k] = fl[k] + min(fl[k], e) + min(max(fl[k] - e, 0) >> 1, R - ft)
492 ]]></artwork>
493 </figure>
494 <t>
495 Here, e is a value that is greater than 0 in case 1, and 0 in case 2.
496 This function is about three times as expensive to evaluate as the
497  high-overhead version, but still an order of magnitude cheaper than a
498  division, since it is composed of very simple operations.
499 </t>
500 <t>
501 In practice it reduces the overhead by about 0.3% of the total bitrate.
502 It also tends to produce R values with a more uniform distribution compared
503  to the high-overhead version, which tends to have peaks in the distribution
504  of R at specific values (see <xref target="SM98"/> for a discussion of this
505  effect).
506 Overall, it makes it more likely that the compression gains from
507  probabilities near one-half are not eliminated by the approximation
508  overhead, increasing the number of symbols that can be usefully modeled.
509 It is an open question whether or not these benefits are worth the increase
510  in computational complexity.
511 </t>
512 </section>
513 </section>
514
515
516 <section title="Context Adaptation">
517 <t>
518 The dyadic adaptation scheme described in&nbsp;<xref target="dyadic_adaptation"/>
519  implements a low-complexity IIR filter for the steady-state case where we only
520  want to adapt the context CDF as fast as the 1/2^r adaptation rate.
521 In many cases, for example when coding symbols at the start of a video frame, only
522  a limited number of symbols have been seen per context.
523 Using this steady-state adaptation scheme risks adapting too slowly and spending
524  too many bits to code symbols with incorrect probability estimates.
525 In other video codecs, this problem is reduced by either implicitly or explicitly
526  allowing for mechanisms to set the initial probability models for a given
527  context.
528 </t>
529 <section title="Implicit Adaptation">
530 <t>
531 One implicit way to use default probabilities is to simply require as a
532  normative part of the decoder that some specific CDFs are used to initialize
533  each context.
534 A representative set of inputs is run through the encoder and a frequency based
535  probability model is computed and reloaded at the start of every frame.
536 This has the advantage of having zero bitstream overhead and is optimal for
537  certain stationary symbols.
538 However for other non-stationary symbols, or highly content dependent contexts
539  where the sample input is not representative, this can be worse than starting
540  with a flat distribution as it now takes even longer to adapt to the
541  steady-state.
542 Moreover the amount of hardware area required to store initial probability
543  tables for each context goes up with the number of contexts in the codec.
544 </t>
545 <t>
546 Another implicit way to deal with poor initial probabilities is through backward
547  adaptation based on the probability estimates from the previous frame.
548 After decoding a frame, the adapted CDFs for each context are simply kept as-is
549  and not reset to their defaults.
550 This has the advantage of having no bitstream overhead, and tracking to certain
551  content types closely as we expect frames with similar content at similar rates,
552  to have well correlated CDFs.
553 However, this only works when we know there will be no bitstream errors due to
554  the transport layer, e.g., TCP or HTTP.
555 In low delay use cases (video on demand, live streaming, video conferencing),
556  implicit backwards adaptation is avoided as it risks desynchronizing the
557  entropy decoder state and permanently losing the video stream.
558 </t>
559 </section>
560 <section title="Explicit Adaptation">
561 <t>
562 For codecs that include the ability to update the probability models in the
563  bitstream, it is possible to explicitly signal a starting CDF.
564 The previously described implicit backwards adaptation is now possible by
565  simply explicitly coding a probability update for each frame.
566 However, the cost of signaling the updated CDF must be overcome by the
567  savings from coding with the updated CDF.
568 Blindly updating all contexts per frame may work at high rates where the size
569  of the CDFs is small relative to the coded symbol data.
570 However at low rates, the benefit of using more accurate CDFs is quickly
571  overcome by the cost of coding them, which increases with the number of
572  contexts.
573 </t>
574 <t>
575 More sophisticated encoders can compute the cost of coding a probability update
576  for a given context, and compare it to the size reduction achieved by coding
577  symbols with this context.
578 Here all symbols for a given frame (or tile) are buffered and not serialized by
579  the entropy coder until the end of the frame (or tile) is reached.
580 Once the end of the entropy segment has been reached, the cost in bits for
581  coding symbols with both the default probabilities and the proposed updated
582  probabilities can be measured and compared.
583 However, note that with the symbols already buffered, rather than consider the
584  context probabilities from the previous frame, a simple frequency based
585  probability model can be computed and measured.
586 Because this probability model is computed based on the symbols we are about
587  to code this technique is called forward adaptation.
588 If the cost in bits to signal and code with this new probability model is less
589  than that of using the default then it is used.
590 This has the advantage of only ever coding a probability update if it is an
591  improvement and producing a bitstream that is robust to errors, but
592  requires an entire entropy segments worth of symbols be cached.
593 </t>
594 </section>
595 <section anchor="early_adaptation" title="Early Adaptation">
596 <t>
597 We would like to take advantage of the low-cost multi-symbol CDF adaptation
598  described in&nbsp;<xref target="dyadic_adaptation"/> without in the broadest set
599  of use cases.
600 This means the initial probability adaptation scheme should support low-delay,
601  error-resilient streams that efficiently implemented in both hardware and
602  software.
603 We propose an early adaptation scheme that supports this goal.
604 </t>
605 <t>
606 At the beginning of a frame (or tile), all CDFs are initialized to a flat
607  distribution.
608 For a given multi-symbol context with M potential symbols, assume that the
609  initial dyadic CDF is initialized so that each symbol has probability 1/M.
610 For the first M coded symbols, the CDF is updated as follows:
611 <figure align="center">
612 <artwork align="center"><![CDATA[
613 a[c,M] = ft/(M + c)
614
615         ( fl[i] - floor((fl[i] - i)*a/ft),          i <= coded symbol
616 fl[i] = <
617         ( fl[i] - floor((fl[i] + M - i - ft)*a/ft), i > coded symbol
618   ]]></artwork>
619 </figure>
620 where c goes from 0 to M-1 and is the running count of the number of symbols
621  coded with this CDF.
622 Note that for a fixed CDF precision (ft is always a power of two) and a
623  maximum number of possible symbols M, the values of a[c,M] can be stored
624  in a M*(M+1)/2 element table, which is 136 entries when M = 16.
625 </t>
626 </section>
627 </section>
628
629 <section anchor="entropy_experiment" title="Simple Experiment">
630 <t>
631 As a simple experiment to validate the non-binary approach, we compared a
632  non-binary arithmetic coder to the VP8 (binary) entropy coder.
633 This was done by instrumenting vp8_treed_read() in libvpx to dump out the
634  symbol decoded and the associated probabilities used to decode it.
635 This data only includes macroblock mode and motion vector information, as the
636  DCT token data is decoded with custom inline functions, and not
637  vp8_treed_read().
638 This data is available at
639  <eref target="https://people.xiph.org/~tterribe/daala/ec_test0/ec_tokens.txt"/>.
640 It includes 1,019,670&nbsp;values encode using 2,125,995&nbsp;binary symbols
641  (or 2.08&nbsp;symbols per value).
642 We expect that with a conscious effort to group symbols during the codec
643  design, this average could easily be increased.
644 </t>
645 <t>
646 We then implemented both the regular VP8 entropy decoder (in plain C, using all
647  of the optimizations available in libvpx at the time) and a multisymbol
648  entropy decoder (also in plain C, using similar optimizations), which encodes
649  each value with a single symbol.
650 For the decoder partition search in the non-binary decoder, we used a simple
651  for loop (O(N) worst-case), even though this could be made constant-time and
652  branchless with a few SIMD instructions such as (on x86) PCMPGTW, PACKUSWB,
653  and PMOVMASKB followed by BSR.
654 The source code for both implementations is available at
655  <eref target="https://people.xiph.org/~tterribe/daala/ec_test0/ec_test.c"/>
656  (compile with -DEC_BINARY for the binary version and -DEC_MULTISYM for the
657  non-binary version).
658 </t>
659 <t>
660 The test simply loads the tokens, and then loops 1024 times encoding them using
661  the probabilities provided, and then decoding them.
662 The loop was added to reduce the impact of the overhead of loading the data,
663  which is implemented very inefficiently.
664 The total runtime on a Core i7 from 2010 is 53.735&nbsp;seconds for the binary
665  version, and 27.937&nbsp;seconds for the non-binary version, or a 1.92x
666  improvement.
667 This is very nearly equal to the number of symbols per value in the binary
668  coder, suggesting that the per-symbol overheads account for the vast majority
669  of the computation time in this implementation.
670 </t>
671 </section>
672
673 </section>
674
675 <section anchor="reversible_integer_transforms"
676  title="Reversible Integer Transforms">
677 <t>
678 Integer transforms in image and video coding date back to at least
679  1969&nbsp;<xref target="PKA69"/>.
680 Although standards such as MPEG2 and MPEG4 Part&nbsp;2 allow some flexibility
681  in the transform implementation, implementations were subject to drift and
682  error accumulation, and encoders had to impose special macroblock refresh
683  requirements to avoid these problems, not always successfully.
684 As transforms in modern codecs only account for on the order of 10% of the
685  total decoder complexity, and, with the use of weighted prediction with gains
686  greater than unity and intra prediction, are far more susceptible to drift and
687  error accumulation, it no longer makes sense to allow a non-exact transform
688  specification.
689 </t>
690 <t>
691 However, it is also possible to make such transforms "reversible", in the sense
692  that applying the inverse transform to the result of the forward transform
693  gives back the original input values, exactly.
694 This gives a lossy codec, which normally quantizes the coefficients before
695  feeding them into the inverse transform, the ability to scale all the way to
696  lossless compression without requiring any new coding tools.
697 This approach has been used successfully by JPEG XR, for
698  example&nbsp;<xref target="TSSRM08"/>.
699 </t>
700 <t>
701 Such reversible transforms can be constructed using "lifting steps", a series
702  of shear operations that can represent any set of plane rotations, and thus
703  any orthogonal transform.
704 This approach dates back to at least 1992&nbsp;<xref target="BE92"/>, which
705  used it to implement a four-point 1-D Discrete Cosine Transform (DCT).
706 Their implementation requires 6&nbsp;multiplications, 10&nbsp;additions,
707  2&nbsp;shifts, and 2&nbsp;negations, and produces output that is a factor of
708  sqrt(2) larger than the orthonormal version of the transform.
709 The expansion of the dynamic range directly translates into more bits to code
710  for lossless compression.
711 Because the least significant bits are usually very nearly random noise, this
712  scaling increases the coding cost by approximately half a bit per sample.
713 </t>
714
715 <section anchor="lifting_steps" title="Lifting Steps">
716 <t>
717 To demonstrate the idea of lifting steps, consider the two-point transform
718 <figure align="center">
719 <artwork align="center"><![CDATA[
720             ___
721 [ y0 ]     / 1  [  1 1 ] [ x0 ]
722 [    ] =  / --- [      ] [    ]
723 [ y1 ]   v   2  [ -1 1 ] [ x1 ]
724 ]]></artwork>
725 </figure>
726 This can be implemented up to scale via
727 <figure align="center">
728 <artwork align="center"><![CDATA[
729 y0 = x0 + x1
730
731 y1 = 2*x1 - y0
732 ]]></artwork>
733 </figure>
734  and reversed via
735 <figure align="center">
736 <artwork align="center"><![CDATA[
737 x1 = (y0 + y1) >> 1
738
739 x0 = y0 - x1
740 ]]></artwork>
741 </figure>
742 </t>
743 <t>
744 Both y0 and y1 are too large by a factor of sqrt(2), however.
745 </t>
746 <t>
747 It is also possible to implement any rotation by an angle t, including the
748  orthonormal scale factor, by decomposing it into three steps:
749 <figure align="center">
750 <artwork align="center"><![CDATA[
751           cos(t) - 1
752 u0 = x0 + ---------- * x1
753             sin(t)
754
755 y1 = x1 + sin(t)*u0
756
757           cos(t) - 1
758 y0 = u0 + ---------- * y1
759             sin(t)
760 ]]></artwork>
761 </figure>
762 By letting t=-pi/4, we get an implementation of the first transform that
763  includes the scaling factor.
764 To get an integer approximation of this transform, we need only replace the
765  transcendental constants by fixed-point approximations:
766 <figure align="center">
767 <artwork align="center"><![CDATA[
768 u0 = x0 + ((27*x1 + 32) >> 6)
769
770 y1 = x1 - ((45*u0 + 32) >> 6)
771
772 y0 = u0 + ((27*y1 + 32) >> 6)
773 ]]></artwork>
774 </figure>
775 This approximation is still perfectly reversible:
776 <figure align="center">
777 <artwork align="center"><![CDATA[
778 u0 = y0 - ((27*y1 + 32) >> 6)
779
780 x1 = y1 + ((45*u0 + 32) >> 6)
781
782 x0 = u0 - ((27*x1 + 32) >> 6)
783 ]]></artwork>
784 </figure>
785 Each of the three steps can be implemented using just two ARM instructions,
786  with constants that have up to 14&nbsp;bits of precision (though using fewer
787  bits allows more efficient hardware implementations, at a small cost in coding
788  gain).
789 However, it is still much more complex than the first approach.
790 </t>
791 <t>
792 We can get a compromise with a slight modification:
793 <figure align="center">
794 <artwork align="center"><![CDATA[
795 y0 = x0 + x1
796
797 y1 = x1 - (y0 >> 1)
798 ]]></artwork>
799 </figure>
800 This still only implements the original orthonormal transform up to scale.
801 The y0 coefficient is too large by a factor of sqrt(2) as before, but y1 is now
802  too small by a factor of sqrt(2).
803 If our goal is simply to (optionally quantize) and code the result, this is
804  good enough.
805 The different scale factors can be incorporated into the quantization matrix in
806  the lossy case, and the total expansion is roughly equivalent to that of the
807  orthonormal transform in the lossless case.
808 Plus, we can perform each step with just one ARM instruction.
809 </t>
810 <t>
811 However, if instead we want to apply additional transformations to the data, or
812  use the result to predict other data, it becomes much more convenient to have
813  uniformly scaled outputs.
814 For a two-point transform, there is little we can do to improve on the
815  three-multiplications approach above.
816 However, for a four-point transform, we can use the last approach and arrange
817  multiple transform stages such that the "too large" and "too small" scaling
818  factors cancel out, producing a result that has the true, uniform, orthonormal
819  scaling.
820 To do this, we need one more tool, which implements the following transform:
821 <figure align="center">
822 <artwork align="center"><![CDATA[
823             ___
824 [ y0 ]     / 1  [ cos(t) -sin(t) ] [ 1  0 ] [ x0 ]
825 [    ] =  / --- [                ] [      ] [    ]
826 [ y1 ]   v   2  [ sin(t)  cos(t) ] [ 0  2 ] [ x1 ]
827 ]]></artwork>
828 </figure>
829 This takes unevenly scaled inputs, rescales them, and then rotates them.
830 Like an ordinary rotation, it can be reduced to three lifting steps:
831 <figure align="center">
832 <artwork align="center"><![CDATA[
833                       _
834           2*cos(t) - v2
835 u0 = x0 + ------------- * x1
836               sin(t)
837             ___
838            / 1
839 y1 = x1 + / --- * sin(t)*u0
840          v   2
841                     _
842           cos(t) - v2
843 y0 = u0 + ----------- * y1
844              sin(t)
845 ]]></artwork>
846 </figure>
847 As before, the transcendental constants may be replaced by fixed-point
848  approximations without harming the reversibility property.
849 </t>
850
851 </section>
852
853 <section anchor="four_point_transform" title="4-Point Transform">
854 <t>
855 Using the tools from the previous section, we can design a reversible integer
856  four-point DCT approximation with uniform, orthonormal scaling.
857 This requires 3&nbsp;multiplies, 9&nbsp;additions, and 2&nbsp;shifts (not
858  counting the shift and rounding offset used in the fixed-point multiplies, as
859  these are built into the multiplier).
860 This is significantly cheaper than the&nbsp;<xref target="BE92"/> approach, and
861  the output scaling is smaller by a factor of sqrt(2), saving half a bit per
862  sample in the lossless case.
863 By comparison, the four-point forward DCT approximation used in VP9, which is
864  not reversible, uses 6&nbsp;multiplies, 6&nbsp;additions, and 2 shifts
865  (counting shifts and rounding offsets which cannot be merged into a single
866  multiply instruction on ARM).
867 Four of its multipliers also require 28-bit accumulators, whereas this proposal
868  can use much smaller multipliers without giving up the reversibility property.
869 The total dynamic range expansion is 1&nbsp;bit: inputs in the range [-256,255)
870  produce transformed values in the range [-512,510).
871 This is the smallest dynamic range expansion possible for any reversible
872  transform constructed from mostly-linear operations.
873 It is possible to make reversible orthogonal transforms with no dynamic range
874  expansion by using "piecewise-linear" rotations&nbsp;<xref target="SLD04"/>,
875  but each step requires a large number of operations in a software
876  implementation.
877 </t>
878
879 <t>
880 Pseudo-code for the forward transform follows:
881 <figure align="left">
882 <artwork align="left"><![CDATA[
883 Input:  x0, x1, x2, x3
884 Output: y0, y1, y2, y3
885 /* Rotate (x3, x0) by -pi/4, asymmetrically scaled output. */
886 t3  = x0 - x3
887 t0  = x0 - (t3 >> 1)
888 /* Rotate (x1, x2) by pi/4, asymmetrically scaled output. */
889 t2  = x1 + x2
890 t2h = t2 >> 1
891 t1  = t2h - x2
892 /* Rotate (t2, t0) by -pi/4, asymmetrically scaled input. */
893 y0  = t0 + t2h
894 y2  = y0 - t2
895 /* Rotate (t3, t1) by 3*pi/8, asymmetrically scaled input. */
896 t3  = t3 - (45*t1 + 32 >> 6)
897 y1  = t1 + (21*t3 + 16 >> 5)
898 y3  = t3 - (71*y1 + 32 >> 6)
899 ]]></artwork>
900 </figure>
901 Even though there are three asymmetrically scaled rotations by pi/4, by careful
902  arrangement we can share one of the shift operations (to help software
903  implementations: shifts by a constant are basically free in hardware).
904 This technique can be used to even greater effect in larger transforms.
905 </t>
906
907 <t>
908 The inverse transform is constructed by simply undoing each step in turn:
909 <figure align="left">
910 <artwork align="left"><![CDATA[
911 Input:  y0, y1, y2, y3
912 Output: x0, x1, x2, x3
913 /* Rotate (y3, y1) by -3*pi/8, asymmetrically scaled output. */
914 t3  = y3 + (71*y1 + 32 >> 6)
915 t1  = y1 - (21*t3 + 16 >> 5)
916 t3  = t3 + (45*t1 + 32 >> 6)
917 /* Rotate (y2, y0) by pi/4, asymmetrically scaled output. */
918 t2  = y0 - y2
919 t2h = t2 >> 1
920 t0  = y0 - t2h
921 /* Rotate (t1, t2) by -pi/4, asymmetrically scaled input. */
922 x2  = t2h - t1
923 x1  = t2 - x2
924 /* Rotate (x3, x0) by pi/4, asymmetrically scaled input. */
925 x0  = t0 - (t3 >> 1)
926 x3  = x0 - t3
927 ]]></artwork>
928 </figure>
929 </t>
930
931 <t>
932 Although the right shifts make this transform non-linear, we can compute
933  "basis functions" for it by sending a vector through it with a single value
934  set to a large constant (256 was used here), and the rest of the values set to
935  zero.
936 The true basis functions for a four-point DCT (up to five digits) are
937 <figure align="left">
938 <artwork align="left"><![CDATA[
939 [ y0 ]   [ 0.50000  0.50000  0.50000  0.50000 ] [ x0 ]
940 [ y1 ] = [ 0.65625  0.26953 -0.26953 -0.65625 ] [ x1 ]
941 [ y2 ]   [ 0.50000 -0.50000 -0.50000  0.50000 ] [ x2 ]
942 [ y3 ]   [ 0.27344 -0.65234  0.65234 -0.27344 ] [ x3 ]
943 ]]></artwork>
944 </figure>
945 The corresponding basis functions for our reversible, integer DCT, computed
946  using the approximation described above, are
947 <figure align="left">
948 <artwork align="left"><![CDATA[
949 [ y0 ]   [ 0.50000  0.50000  0.50000  0.50000 ] [ x0 ]
950 [ y1 ] = [ 0.65328  0.27060 -0.27060 -0.65328 ] [ x1 ]
951 [ y2 ]   [ 0.50000 -0.50000 -0.50000  0.50000 ] [ x2 ]
952 [ y3 ]   [ 0.27060 -0.65328  0.65328 -0.27060 ] [ x3 ]
953 ]]></artwork>
954 </figure>
955 The mean squared error (MSE) of the output, compared to a true DCT, can be
956  computed with some assumptions about the input signal.
957 Let G be the true DCT basis and G' be the basis for our integer approximation
958  (computed as described above).
959 Then the error in the transformed results is
960 <figure align="left">
961 <artwork align="left"><![CDATA[
962 e = G.x - G'.x = (G - G').x = D.x
963 ]]></artwork>
964 </figure>
965  where D&nbsp;=&nbsp;(G&nbsp;-&nbsp;G')&nbsp;.
966 The MSE is then&nbsp;<xref target="Que98"/>
967 <figure align="left">
968 <artwork align="left"><![CDATA[
969 1              1
970 - * E[e^T.e] = - * E[x^T.D^T.D.x]
971 N              N
972
973                1
974              = - * E[tr(D.x.x^T.D^T)]
975                N
976
977                1
978              = - * E[tr(D.Rxx.D^T)]
979                N
980 ]]></artwork>
981 </figure>
982  where Rxx is the autocorrelation matrix of the input signal.
983 Assuming the input is a zero-mean, first-order autoregressive (AR(1)) process
984  gives an autocorrelation matrix of
985 <figure align="left">
986 <artwork align="left"><![CDATA[
987               |i - j|
988 Rxx[i,j] = rho
989 ]]></artwork>
990 </figure>
991  for some correlation coefficient rho.
992 A value of rho&nbsp;=&nbsp;0.95 is typical for image compression applications.
993 Smaller values are more normal for motion-compensated frame differences, but
994  this makes surprisingly little difference in transform design.
995 Using the above procedure, the theoretical MSE of this approximation is
996  1.230E-6, which is below the level of the truncation error introduced by the
997  right shift operations.
998 This suggests the dynamic range of the input would have to be more than
999  20&nbsp;bits before it became worthwhile to increase the precision of the
1000  constants used in the multiplications to improve accuracy, though it may be
1001  worth using more precision to reduce bias.
1002 </t>
1003
1004 </section>
1005
1006 <section anchor="larger_transforms" title="Larger Transforms">
1007 <t>
1008 The same techniques can be applied to construct a reversible eight-point DCT
1009  approximation with uniform, orthonormal scaling using 15&nbsp;multiplies,
1010  31&nbsp;additions, and 5&nbsp;shifts.
1011 It is possible to reduce this to 11&nbsp;multiplies and 29&nbsp;additions,
1012  which is the minimum number of multiplies possible for an eight-point DCT with
1013  uniform scaling&nbsp;<xref target="LLM89"/>, by introducing a scaling factor
1014  of sqrt(2), but this harms lossless performance.
1015 The dynamic range expansion is 1.5&nbsp;bits (again the smallest possible), and
1016  the MSE is 1.592E-06.
1017 By comparison, the eight-point transform in VP9 uses 12&nbsp;multiplications,
1018  32&nbsp;additions, and 6 shifts.
1019 </t>
1020 <t>
1021 Similarly, we have constructed a reversible sixteen-point DCT approximation
1022  with uniform, orthonormal scaling using 33&nbsp;multiplies, 83&nbsp;additions,
1023  and 16&nbsp;shifts.
1024 This is just 2&nbsp;multiplies and 2&nbsp;additions more than the
1025  (non-reversible, non-integer, but uniformly scaled) factorization
1026  in&nbsp;<xref target="LLM89"/>.
1027 By comparison, the sixteen-point transform in VP9 uses 44&nbsp;multiplies,
1028  88&nbsp;additions, and 18&nbsp;shifts.
1029 The dynamic range expansion is only 2&nbsp;bits (again the smallest possible),
1030  and the MSE is 1.495E-5.
1031 </t>
1032 <t>
1033 We also have a reversible 32-point DCT approximation with uniform,
1034  orthonormal scaling using 87&nbsp;multiplies, 215&nbsp;additions, and
1035  38&nbsp;shifts.
1036 By comparison, the 32-point transform in VP9 uses 116&nbsp;multiplies,
1037  194&nbsp;additions, and 66&nbsp;shifts.
1038 Our dynamic range expansion is still the minimal 2.5&nbsp;bits, and the MSE is
1039  8.006E-05
1040 </t>
1041 <t>
1042 Code for all of these transforms is available in the development repository
1043  listed in&nbsp;<xref target="development_repository"/>.
1044 </t>
1045 </section>
1046
1047 <section anchor="hadamard_transforms" title="Walsh-Hadamard Transforms">
1048 <t>
1049 These techniques can also be applied to constructing Walsh-Hadamard
1050  Transforms, another useful transform family that is cheaper to implement than
1051  the DCT (since it requires no multiplications at all).
1052 The WHT has many applications as a cheap way to approximately change the time
1053  and frequency resolution of a set of data (either individual bands, as in the
1054  Opus audio codec, or whole blocks).
1055 VP9 uses it as a reversible transform with uniform, orthonormal scaling for
1056  lossless coding in place of its DCT, which does not have these properties.
1057 </t>
1058
1059 <t>
1060 Applying a 2x2 WHT to a block of 2x2 inputs involves running a 2-point WHT on
1061  the rows, and then another 2-point WHT on the columns.
1062 The basis functions for the 2-point WHT are, up to scaling, [1,&nbsp;1] and
1063  [1,&nbsp;-1].
1064 The four variations of a two-step lifer given in
1065  <xref target="lifting_steps"/> are exactly the lifting steps needed to
1066  implement a 2x2 WHT: two stages that produce asymmetrically scaled outputs
1067  followed by two stages that consume asymmetrically scaled inputs.
1068 <figure align="left">
1069 <artwork align="left"><![CDATA[
1070 Input:  x00, x01, x10, x11
1071 Output: y00, y01, y10, y11
1072 /* Transform rows */
1073 t1 = x00 - x01
1074 t0 = x00 - (t1 >> 1) /* == (x00 + x01)/2 */
1075 t2 = x10 + x11
1076 t3 = (t2 >> 1) - x11 /* == (x10 - x11)/2 */
1077 /* Transform columns */
1078 y00 = t0 + (t2 >> 1) /* == (x00 + x01 + x10 + x11)/2 */
1079 y10 = y00 - t2       /* == (x00 + x01 - x10 - x11)/2 */
1080 y11 = (t1 >> 1) - t3 /* == (x00 - x01 - x10 + x11)/2 */
1081 y01 = t1 - y11       /* == (x00 - x01 + x10 - x11)/2 */
1082 ]]></artwork>
1083 </figure>
1084 </t>
1085
1086 <t>
1087 By simply re-ordering the operations, we can see that there are two shifts that
1088  may be shared between the two stages:
1089 <figure align="left">
1090 <artwork align="left"><![CDATA[
1091 Input:  x00, x01, x10, x11
1092 Output: y00, y01, y10, y11
1093 t1 = x00 - x01
1094 t2 = x10 + x11
1095 t0 = x00 - (t1 >> 1) /* == (x00 + x01)/2 */
1096 y00 = t0 + (t2 >> 1) /* == (x00 + x01 + x10 + x11)/2 */
1097 t3 = (t2 >> 1) - x11 /* == (x10 - x11)/2 */
1098 y11 = (t1 >> 1) - t3 /* == (x00 - x01 - x10 + x11)/2 */
1099 y10 = y00 - t2       /* == (x00 + x01 - x10 - x11)/2 */
1100 y01 = t1 - y11       /* == (x00 - x01 + x10 - x11)/2 */
1101 ]]></artwork>
1102 </figure>
1103 </t>
1104
1105 <t>
1106 By eliminating the double-negation of x11 and re-ordering the additions to it,
1107  we can see even more operations in common:
1108 <figure align="left">
1109 <artwork align="left"><![CDATA[
1110 Input:  x00, x01, x10, x11
1111 Output: y00, y01, y10, y11
1112 t1 = x00 - x01
1113 t2 = x10 + x11
1114 t0 = x00 - (t1 >> 1) /* == (x00 + x01)/2 */
1115 y00 = t0 + (t2 >> 1) /* == (x00 + x01 + x10 + x11)/2 */
1116 t3 = x11 + (t1 >> 1) /* == x11 + (x00 - x01)/2 */
1117 y11 = t3 - (t2 >> 1) /* == (x00 - x01 - x10 + x11)/2 */
1118 y10 = y00 - t2       /* == (x00 + x01 - x10 - x11)/2 */
1119 y01 = t1 - y11       /* == (x00 - x01 + x10 - x11)/2 */
1120 ]]></artwork>
1121 </figure>
1122 </t>
1123
1124 <t>
1125 Simplifying further, the whole transform may be computed with just
1126  7&nbsp;additions and 1&nbsp;shift:
1127 <figure align="left">
1128 <artwork align="left"><![CDATA[
1129 Input:  x00, x01, x10, x11
1130 Output: y00, y01, y10, y11
1131 t1 = x00 - x01
1132 t2 = x10 + x11
1133 t4 = (t2 - t1) >> 1 /* == (-x00 + x01 + x10 + x11)/2 */
1134 y00 = x00 + t4      /* ==  (x00 + x01 + x10 + x11)/2 */
1135 y11 = x11 - t4      /* ==  (x00 - x01 - x10 + x11)/2 */
1136 y10 = y00 - t2      /* ==  (x00 + x01 - x10 - x11)/2 */
1137 y01 = t1 - y11      /* ==  (x00 - x01 + x10 - x11)/2 */
1138 ]]></artwork>
1139 </figure>
1140 </t>
1141
1142 <t>
1143 This is a significant savings over other approaches described in the
1144  literature, which require 8&nbsp;additions, 2&nbsp;shifts, and
1145  1&nbsp;negation&nbsp;<xref target="FOIK99"/> (37.5%&nbsp;more operations), or
1146  10&nbsp;additions, 1&nbsp;shift, and
1147  2&nbsp;negations&nbsp;<xref target="TSSRM08"/> (62.5%&nbsp;more operations).
1148 The same operations can be applied to compute a 4-point WHT in one dimension.
1149 This implementation is used in this way in VP9's lossless mode.
1150 Since larger WHTs may be trivially factored into multiple smaller WHTs, the
1151  same approach can implement a reversible, orthonormally scaled WHT of any size
1152  (2**N)x(2**M), so long as (N&nbsp;+&nbsp;M) is even.
1153 </t>
1154
1155 </section>
1156
1157 </section>
1158
1159 <section anchor="development_repository" title="Development Repository">
1160 <t>
1161 The tools presented here were developed as part of Xiph.Org's Daala project.
1162 They are available, along with many others in greater and lesser states of
1163  maturity, in the Daala git repository at
1164  <eref target="https://git.xiph.org/daala.git"/>.
1165 See <eref target="https://xiph.org/daala/"/> for more information.
1166 </t>
1167 </section>
1168
1169 <section title="IANA Considerations">
1170 <t>
1171 This document has no actions for IANA.
1172 </t>
1173 </section>
1174
1175 <section anchor="Acknowledgments" title="Acknowledgments">
1176 <t>
1177 Thanks to Nathan Egge, Gregory Maxwell, and Jean-Marc Valin for their
1178  assistance in the implementation and experimentation, and in preparing this
1179  draft.
1180 </t>
1181 </section>
1182
1183 </middle>
1184 <back>
1185 <!--references title="Normative References">
1186
1187 <?rfc include="http://xml.resource.org/public/rfc/bibxml/reference.RFC.2119.xml"?>
1188
1189 </references-->
1190
1191 <references title="Informative References">
1192
1193 <?rfc include="http://xml.resource.org/public/rfc/bibxml/reference.RFC.6386.xml"?>
1194 <?rfc include="http://xml.resource.org/public/rfc/bibxml/reference.RFC.6716.xml"?>
1195
1196 <reference anchor="BE92">
1197 <front>
1198 <title>New Networks for Perfect Inversion and Perfect Reconstruction</title>
1199 <author initials="F.A.M.L." surname="Bruekers"
1200  fullname="Fons A.M.L. Bruekers"/>
1201 <author initials="A.W.M." surname="van den Enden"
1202  fullname="Ad W.M. van den Enden"/>
1203 <date month="January" year="1992"/>
1204 </front>
1205 <seriesInfo name="IEEE Journal on Selected Areas in Communication"
1206  value="10(1):129--137"/>
1207 </reference>
1208
1209 <reference anchor="FOIK99">
1210 <front>
1211 <title>Lossless 8-point Fast Discrete Cosine Transform Using Lossless
1212  Hadamard Transform</title>
1213 <author initials="S." surname="Fukuma" fullname="Shinji Fukuma"/>
1214 <author initials="K." surname="Oyama" fullname="Koichi Oyama"/>
1215 <author initials="M." surname="Iwahashi" fullname="Masahiro Iwahashi"/>
1216 <author initials="N." surname="Kambayashi" fullname="Noriyoshi Kambayashi"/>
1217 <date month="October" year="1999"/>
1218 </front>
1219 <seriesInfo name="Technical Report"
1220  value="The Institute of Electronics, Information, and Communication Engineers
1221  of Japan"/>
1222 </reference>
1223
1224 <reference anchor="LLM89">
1225 <front>
1226 <title>Practical Fast 1-D DCT Algorithms with 11 Multiplications</title>
1227 <author initials="C." surname="Loeffler" fullname="Christoph Loeffler"/>
1228 <author initials="A." surname="Ligtenberg" fullname="Adriaan Ligtenberg"/>
1229 <author initials="G.S." surname="Moschytz" fullname="George S. Moschytz"/>
1230 <date month="May" year="1989"/>
1231 </front>
1232 <seriesInfo name="Proc. Acoustics, Speech, and Signal Processing (ICASSP'89)"
1233  value="vol. 2, pp. 988--991"/>
1234 </reference>
1235
1236 <reference anchor="Pas76">
1237 <front>
1238 <title>Source Coding Algorithms for Fast Data Compression</title>
1239 <author initials="R.C." surname="Pasco" fullname="Richard C. Pasco"/>
1240 <date month="May" year="1976"/>
1241 </front>
1242 <seriesInfo name="Ph.D. Thesis"
1243  value="Dept. of Electrical Engineering, Stanford University"/>
1244 </reference>
1245
1246 <reference anchor="PKA69">
1247 <front>
1248 <title>Hadamard Transform Image Coding</title>
1249 <author initials="W.K." surname="Pratt" fullname="W.K. Pratt"/>
1250 <author initials="J." surname="Kane" fullname="J. Kane"/>
1251 <author initials="H.C." surname="Andrews" fullname="H.C. Andrews"/>
1252 <date month="Jan" year="1969"/>
1253 </front>
1254 <seriesInfo name="Proc. IEEE" value="57(1):58--68"/>
1255 </reference>
1256
1257 <reference anchor="Que98">
1258 <front>
1259 <title>On Unitary Transform Approximations</title>
1260 <author initials="R.L." surname="de Queiroz" fullname="Ricardo L. de Queiroz"/>
1261 <date month="Feb" year="1998"/>
1262 </front>
1263 <seriesInfo name="IEEE Signal Processing Letters" value="5(2):46--47"/>
1264 </reference>
1265
1266 <reference anchor="SLD04">
1267 <front>
1268 <title>An Improved N-Bit to N-Bit Reversible Haar-Like Transform</title>
1269 <author initials="J.G." surname="Senecal" fullname="Joshua G. Senecal"/>
1270 <author initials="P." surname="Lindstrom" fullname="Peter Lindstrom"/>
1271 <author initials="M.A." surname="Duchaineau" fullname="Mark A. Duchaineau"/>
1272 <date month="October" year="2004"/>
1273 </front>
1274 <seriesInfo
1275  name="Proc. of the 12th Pacific Conference on Computer Graphics and Applications (PG'04)"
1276  value="pp. 371--380"/>
1277 </reference>
1278
1279 <reference anchor="SM98">
1280 <front>
1281 <title>Piecewise Integer Mapping for Arithmetic Coding</title>
1282 <author initials="L." surname="Stuiver" fullname="Lang Stuiver"/>
1283 <author initials="A." surname="Moffat" fullname="Alistair Moffat"/>
1284 <date month="March/April" year="1998"/>
1285 </front>
1286 <seriesInfo name="Proc. of the 17th IEEE Data Compression Conference (DCC'98)"
1287  value="pp. 1--10"/>
1288 </reference>
1289
1290 <reference anchor="TSSRM08">
1291 <front>
1292 <title>Low-complexity Hierarchical Lapped Transform for Lossy-to-Lossless
1293  Image Coding in JPEG XR/HD Photo</title>
1294 <author initials="C." surname="Tu" fullname="Chengjie Tu"/>
1295 <author initials="S." surname="Srinivasan" fullname="Sridhar Srinivasan"/>
1296 <author initials="G.J." surname="Sullivan" fullname="Gary J. Sullivan"/>
1297 <author initials="S." surname="Regunathan" fullname="Shankar Regunathan"/>
1298 <author initials="H.S." surname="Malvar" fullname="Henrique S. Malvar"/>
1299 <date month="August" year="2008"/>
1300 </front>
1301 <seriesInfo name="Applications of Digital Image Processing XXXI"
1302  value="vol 7073"/>
1303 </reference>
1304
1305 </references>
1306
1307 </back>
1308 </rfc>