Off[General::spell1]

<< Graphics`Graphics`

First we set up a few contrast sensitivity functions using the data from chapter 5 of Nadenau's thesis. There is an error in Table 5.3 on page 83: the roles of a1 and a2 are reversed. Simply swapping the data for these two parameters yields correct results, i.e. the graph below matches the one in Figure 5.9 on page 85.

sl[f_, a1_, a2_, b1_, b2_, c1_, c2_, s_] := 
      s*(a1*(f^2)*Exp[b1*(f^c1)] + a2*Exp[b2*(f^c2)])

sc[f_, a1_, b1_, c1_, s_] := s*(a1*Exp[b1*(f^c1)])

ycbcrsl[f_] := 
      sl[f, 0.9973, 0.221, -0.9699, -0.800, 0.7578, 1.999, 267.5]

ycbcrsrg[f_] := sc[f, 1.000, -0.1521, 0.893, 45.44]

ycbcrsby[f_] := sc[f, 1.00, -0.2041, 0.900, 28.92]

LogLogPlot[{ycbcrsl[f], ycbcrsrg[f], ycbcrsby[f]}, {f, .5, 20}, 
      PlotRange -> {{0.5, 20}, {10, 300}}, AspectRatio -> 2, 
      PlotStyle -> {{RGBColor[0, 0, 0]}, {RGBColor[1, 0, 0]}, {
            RGBColor[0, 0, 1]}}];

Now we set up variations of the CSF functions for different contrast levels.

ycbcrstsl[f_, c_] := 1/(c^(1/.5) + 1/ycbcrsl[f])

ycbcrstsrg[f_, c_] := 1/(c^(1/.5) + 1/ycbcrsrg[f])

ycbcrstsby[f_, c_] := 1/(c^(1/.5) + 1/ycbcrsby[f])

LogLogPlot[{ycbcrstsl[f, .1], ycbcrstsl[f, .2], ycbcrstsl[f, .3], 
        ycbcrstsl[f, .4], ycbcrstsl[f, .5], ycbcrstsl[f, .6], 
        ycbcrstsl[f, .7], ycbcrstsl[f, .8], ycbcrstsl[f, .9], 
        ycbcrstsl[f, 1.0]}, {f, .5, 30}, PlotRange -> {{.5, 30}, {.5, 75}}, 
      AspectRatio -> 2, PlotStyle -> {{RGBColor[0, 0, 0]}}];

LogLogPlot[{ycbcrstsrg[f, .1], ycbcrstsrg[f, .2], ycbcrstsrg[f, .3], 
        ycbcrstsrg[f, .4], ycbcrstsrg[f, .5], ycbcrstsrg[f, .6], 
        ycbcrstsrg[f, .7], ycbcrstsrg[f, .8], ycbcrstsrg[f, .9], 
        ycbcrstsrg[f, 1.0]}, {f, .5, 30}, PlotRange -> {{.5, 30}, {.5, 30}}, 
      AspectRatio -> 2, PlotStyle -> {{RGBColor[1, 0, 0]}}];

LogLogPlot[{ycbcrstsby[f, .1], ycbcrstsby[f, .2], ycbcrstsby[f, .3], 
        ycbcrstsby[f, .4], ycbcrstsby[f, .5], ycbcrstsby[f, .6], 
        ycbcrstsby[f, .7], ycbcrstsby[f, .8], ycbcrstsby[f, .9], 
        ycbcrstsby[f, 1.0]}, {f, .5, 30}, PlotRange -> {{.5, 30}, {.5, 21}}, 
      AspectRatio -> 2, PlotStyle -> {{RGBColor[0, 0, 1]}}];

In order to ensure that increasing the viewing distance does not expose artifacts, we need to boost the luma CSF curve to its maximum value at all frequencies less than the one where it achieves that maximum. The chroma CSF curves have their maximum at frequency 0.

ycbcrslpeakf = 
      ReplaceAll[f, FindMinimum[-ycbcrsl[f], {f, 5}][[2]]]

Now we normalize all the curves so their maximums are 1.

ycbcrsln[f_] := 
      If[f < ycbcrslpeakf, 1, 
        sl[f, 0.9973, 0.221, -0.9699, -0.800, 0.7578, 1.999, 1.0]/
          sl[ycbcrslpeakf, 0.9973, 0.221, -0.9699, -0.800, 0.7578, 
            1.999, 1.0]]

ycbcrsrgn[f_] := 
      sc[f, 1.000, -0.1521, 0.893, 45.44]/
        sc[0, 1.000, -0.1521, 0.893, 45.44]

ycbcrsbyn[f_] := 
      sc[f, 1.00, -0.2041, 0.900, 28.92]/
        sc[0, 1.00, -0.2041, 0.900, 28.92]

Make some nice compiled versions of these functions, or stuff below will become unbearably slow.

ycbcrslnc = Compile[{f}, Evaluate[ycbcrsln[f]]];

ycbcrsrgnc = Compile[{f}, Evaluate[ycbcrsrgn[f]]];

ycbcrsbync = Compile[{f}, Evaluate[ycbcrsbyn[f]]];

ycbcrstsln[f_, c_] := 
      (c^(1/.5) + 1/ycbcrsln[ycbcrslpeakf])/
        (c^(1/.5) + 1/ycbcrsln[f])

ycbcrstsrgn[f_, c_] := 
      (c^(1/.5) + 1/ycbcrsrgn[0])/
        (c^(1/.5) + 1/ycbcrsrgn[f])

ycbcrstsbyn[f_, c_] := 
      (c^(1/.5) + 1/ycbcrsbyn[0])/
        (c^(1/.5) + 1/ycbcrsbyn[f])

ycbcrstslnc = Compile[{f, c}, Evaluate[ycbcrstsln[f, c]]];

ycbcrstsrgnc = Compile[{f, c}, Evaluate[ycbcrstsrgn[f, c]]];

ycbcrstsbync = Compile[{f, c}, Evaluate[ycbcrstsbyn[f, c]]];

fmax[v_, r_] := v*r*(Tan[π/360]*10000/254)

fmaxcrtl := fmax[.50, 90]

fmaxcrtc := fmaxcrtl/2

Plot[{ycbcrstslnc[f, 0], ycbcrstslnc[f, 1], ycbcrstslnc[f, 2], 
        ycbcrstslnc[f, 4], ycbcrstslnc[f, 8], ycbcrstslnc[f, 16], 
        ycbcrstslnc[f, 32], ycbcrstslnc[f, 64]}, {f, 0, fmaxcrtl}, 
      PlotRange -> {{0, fmaxcrtl}, {0, 1}}, AspectRatio -> 2, 
      PlotStyle -> {{RGBColor[0, 0, 0]}}];

Plot[{ycbcrstsrgnc[f, 0], ycbcrstsrgnc[f, 1], ycbcrstsrgnc[f, 2], 
        ycbcrstsrgnc[f, 4], ycbcrstsrgnc[f, 8], ycbcrstsrgnc[f, 16], 
        ycbcrstsrgnc[f, 32], ycbcrstsrgnc[f, 64]}, {f, 0, fmaxcrtc}, 
      PlotRange -> {{0, fmaxcrtc}, {0, 1}}, AspectRatio -> 2, 
      PlotStyle -> {{RGBColor[1, 0, 0]}}];

Plot[{ycbcrstsbync[f, 0], ycbcrstsbync[f, 1], ycbcrstsbync[f, 2], 
        ycbcrstsbync[f, 4], ycbcrstsbync[f, 8], ycbcrstsbync[f, 16], 
        ycbcrstsbync[f, 32], ycbcrstsbync[f, 64]}, {f, 0, fmaxcrtc}, 
      PlotRange -> {{0, fmaxcrtc}, {0, 1}}, AspectRatio -> 2, 
      PlotStyle -> {{RGBColor[0, 0, 1]}}];

<< Graphics`MultipleListPlot`

trapezoidal[l_, h_] := 
      (h/(Length[l] - 1))*
        ((1/2)*(l[[1]] + l[[Length[l]]]) + 
            (Plus@@Take[l, {2, Length[l] - 1}]))

simpson[l_, h_] := 
      (h/(3*Length[l] - 3))*
        (Plus@@
            MapIndexed[#1*
                  (If[EvenQ[#], 4, If[# > 1 && # < Length[l], 2, 1]]&
                        [#2[[1]]])&, l])

mag2[v_] := Plus@@(#*#&/@v)

mag[v_] := Sqrt[mag2[v]]

fwtable[a_, n_] := 
      Sqrt[Length[a]]*
        Transpose[
          #/mag2[#]&/@
            Transpose[
              Sqrt[2*n]*
                (Abs[
                        Drop[Fourier[Join[#, Table[0, {2*n - Length[#]}]]], 
                          -n + 1]]&/@a)]]

We define Theora's DCT in terms of its iDCT implementation, since that is standardized while the fDCT implementation is not. The values below were obtained by pretending the integer transform was linear, and replacing all integer operations with floating point ones in order to compute the basis functions.

theoraidct8
       := {{0.176776885986328125, 0.245197296142578125, 0.230968475341796875, 
          0.207866668701171875, 0.176776885986328125, 0.13889312744140625, 
          0.095672607421875, 0.048770904541015625}, {0.176776885986328125, 
          0.20786713238339871168136596679688, 0.095672607421875, 
          -0.04877171135740354657173156738281, -0.176776885986328125, 
          -0.24519646778935566544532775878906, -0.230968475341796875, 
          -0.13889458333142101764678955078125}, {0.176776885986328125, 
          0.13889458333142101764678955078125, -0.095672607421875, 
          -0.24519646778935566544532775878906, -0.176776885986328125, 
          0.04877171135740354657173156738281, 0.230968475341796875, 
          0.20786713238339871168136596679688}, {0.176776885986328125, 
          0.048770904541015625, -0.230968475341796875, 
          -0.13889312744140625, 0.176776885986328125, 
          0.207866668701171875, -0.095672607421875, 
          -0.245197296142578125}, {0.176776885986328125, 
          -0.048770904541015625, -0.230968475341796875, 
          0.13889312744140625, 0.176776885986328125, 
          -0.207866668701171875, -0.095672607421875, 
          0.245197296142578125}, {0.176776885986328125, 
          -0.13889458333142101764678955078125, -0.095672607421875, 
          0.24519646778935566544532775878906, -0.176776885986328125, 
          -0.04877171135740354657173156738281, 0.230968475341796875, 
          -0.20786713238339871168136596679688}, {0.176776885986328125, 
          -0.20786713238339871168136596679688, 0.095672607421875, 
          0.04877171135740354657173156738281, -0.176776885986328125, 
          0.24519646778935566544532775878906, -0.230968475341796875, 
          0.13889458333142101764678955078125}, {0.176776885986328125, 
          -0.245197296142578125, 0.230968475341796875, 
          -0.207866668701171875, 0.176776885986328125, 
          -0.13889312744140625, 0.095672607421875, 
          -0.048770904541015625}}

theoradct8 := Inverse[theoraidct8]

Show[
      GraphicsArray[
        Table[BarChart[theoradct8[[i + j*2]], Axes -> {False, True}, 
            BarGroupSpacing -> 0, BarStyle -> {GrayLevel[0]}, 
            PlotRange -> {{0, 8}, {-1, 1}}, 
            DisplayFunction -> Identity], {j, 0, 3}, {i, 2}]]];

Compute a table of weighted frequency responses for each coefficient.

theoradct8fwtable1024 = fwtable[theoradct8, 1024];

ListPlot[#, PlotJoined -> True, Axes -> {False, True}]&/@
      theoradct8fwtable1024;

The computation of the CSF filters themselves has been moved to tools/csfgen.c
The version that used to reside in this notebook was fundamentally flawed, and the version implemented in that file has not been ported back into Mathematica. Correspondingly, the code to dump the tables to a C file has been disabled.

Finally, write these to a C file for cut-and-pasting purposes.