Notebook accompanying the paper
How were the MP3 filter coeffients produced?

by Tom H. Koornwinder, version of 5 February 2007, http://www.science.uva.nl/~thk/art/misc/

Korteweg-de Vries Institute, University of Amsterdam, thk@science.uva.nl

The analysis filter coefficients C[i] (i = 1,…,255) taken from Annex C, Table C.1 in:
Information technology --- Coding of moving pictures and associated audio for digital storage media at up to about 1,5 Mbit/s --- Part 3: Audio,
ISO/IEC 11172-3:1993, International Organization for Standardization, 1993;
http://www.iso.org/iso/en/prods-services/popstds/mpeg.html  or
http://www.cise.ufl.edu/~jfd/mp3links.html.
The coefficients are here written as c[[i]].

c = {-0.000000477, -0.000000477, -0.000000477, -0.000000477, -0.000000477, -0.000000477, -0.000000954, -0.000000954, -0.000000954, -0.000000954, -0.000001431, -0.000001431, -0.000001907, -0.000001907, -0.000002384, -0.000002384, -0.000002861, -0.000003338, -0.000003338, -0.000003815, -0.000004292, -0.000004768, -0.000005245, -0.000006199, -0.000006676, -0.000007629, -0.000008106, -0.000009060, -0.000010014, -0.000011444, -0.000012398, -0.000013828, -0.000014782, -0.000016689, -0.000018120, -0.000019550, -0.000021458, -0.000023365, -0.000025272, -0.000027657, -0.000030041, -0.000032425, -0.000034809, -0.000037670, -0.000040531, -0.000043392, -0.000046253, -0.000049591, -0.000052929, -0.000055790, -0.000059605, -0.000062943, -0.000066280, -0.000070095, -0.000073433, -0.000076771, -0.000080585, -0.000083923, -0.000087261, -0.000090599, -0.000093460, -0.000096321, -0.000099182, 
          0.000101566, 0.000103951, 0.000105858, 0.000107288, 0.000108242, 
          0.000108719, 0.000108719, 0.000108242, 0.000106812, 0.000105381, 
          0.000102520, 0.000099182, 0.000095367, 0.000090122, 0.000084400, 
          0.000077724, 0.000069618, 0.000060558, 0.000050545, 0.000039577, 
          0.000027180, 
          0.000013828, -0.000000954, -0.000017166, -0.000034332, -0.000052929, -0.000072956, -0.000093937, -0.000116348, -0.000140190, -0.000165462, -0.000191212, -0.000218868, -0.000247478, -0.000277042, -0.000307560, -0.000339031, -0.000371456, -0.000404358, -0.000438213, -0.000472546, -0.000507355, -0.000542164, -0.000576973, -0.000611782, -0.000646591, -0.000680923, -0.000714302, -0.000747204, -0.000779152, -0.000809669, -0.000838757, -0.000866413, -0.000891685, -0.000915051, -0.000935555, -0.000954151, -0.000968933, -0.000980854, -0.000989437, -0.000994205, -0.000995159, -0.000991821, -0.000983715, 0.000971317, 
          0.000953674, 0.000930786, 0.000902653, 0.000868797, 0.000829220, 
          0.000783920, 0.000731945, 0.000674248, 0.000610352, 0.000539303, 
          0.000462532, 0.000378609, 0.000288486, 0.000191689, 
          0.000088215, -0.000021458, -0.000137329, -0.000259876, -0.000388145, -0.000522137, -0.000661850, -0.000806808, -0.000956535, -0.001111031, -0.001269817, -0.001432419, -0.001597881, -0.001766682, -0.001937389, -0.002110004, -0.002283096, -0.002457142, -0.002630711, -0.002803326, -0.002974033, -0.003141880, -0.003306866, -0.003467083, -0.003622532, -0.003771782, -0.003914356, -0.004048824, -0.004174709, -0.004290581, -0.004395962, -0.004489899, -0.004570484, -0.004638195, -0.004691124, -0.004728317, -0.004748821, -0.004752159, -0.004737377, -0.004703045, -0.004649162, -0.004573822, -0.004477024, -0.004357815, -0.004215240, -0.004049301, -0.003858566, -0.003643036, -0.003401756, 0.003134727, 0.002841473, 0.002521515, 0.002174854, 
          0.001800537, 0.001399517, 0.000971317, 0.000515938, 
          0.000033379, -0.000475883, -0.001011848, -0.001573563, -0.002161503, -0.002774239, -0.003411293, -0.004072189, -0.004756451, -0.005462170, -0.006189346, -0.006937027, -0.007703304, -0.008487225, -0.009287834, -0.010103703, -0.010933399, -0.011775017, -0.012627602, -0.013489246, -0.014358521, -0.015233517, -0.016112804, -0.016994476, -0.017876148, -0.018756866, -0.019634247, -0.020506859, -0.021372318, -0.022228718, -0.023074150, -0.023907185, -0.024725437, -0.025527000, -0.026310921, -0.027073860, -0.027815342, -0.028532982, -0.029224873, -0.029890060, -0.030526638, -0.031132698, -0.031706810, -0.032248020, -0.032754898, -0.033225536, -0.033659935, -0.034055710, -0.034412861, -0.034730434, -0.035007000, -0.035242081, -0.035435200, -0.035586357, -0.035694122, -0.035758972};

From the full list in the cited reference we can observe a symmetry, by which the list can be completed as follows:

c = Flatten[{c, 0.035780907, Reverse[-c]}];

From this list the real (smooth) filter coefficients (Subscript[h, n]) (here written as h[[n]]) can be produced:

For[k = 5, k ≤ 7, k++, 
        c[[64  k]] = -c[[64 k]]];
h = 
        Table[(-1)^(Quotient[n, 64]) c[[n]]/2, {n, 1, 
            511}];
ListPlot[h, PlotRange -> {All, All}];

Then the  (Subscript[h, n])are symmetric:

(1a)     (Subscript[h, n])= (Subscript[h, 512 - n])

The cosine polynomial  φ(ξ) is essentially the Fourier transform of the (Subscript[h, i]):

φ[ξ_] := 
      h[[256]] + 
        2  Sum[h[[256 - n]] Cos[n ξ], {n, 1, 255}]

This is 2π-periodic and even:

(1b)     φ(ξ) = φ(ξ+2π) ,     φ(ξ)  = φ(– ξ)

The graph of φ[ξ] on the interval [0, π / 32] is shown below. The graph in red is the ideal Fourier transform, which is not attainable with finitely many filter coefficients.

Plot[{φ[ξ], UnitStep[π/64 - ξ]}, {ξ, 
          0, π/32}, PlotStyle -> {Black, Red}, 
        Ticks -> {{π/128, π/64, 3  π/128, π/32}, 
            Automatic}];

The graph of φ(ξ) on the interval [π/32, π] is shown below. There φ(ξ) almost vanishes.

Plot[φ[ξ], {ξ, π/32, π}, 
        Ticks -> {{π/32, π}, Automatic}];

For optimal filter properties we should have the normalization

(2a)    (1/(2  π))(∫_(-π)^π) ((φ(ξ))^2) dξ ≈ (1/64)

or equivalently,

(2b)     \(TraditionalForm\`\(\(\ \)\(∑\)\)\+\(\(\ \ \)\(n = 1\)\)\%\(\(\ \ \)\(511\)\)\)  (Subscript[h, n]^(2 ))= h·h ≈ (1/64)

This is indeed satisfied:

h . h - (1/64)

A further requirement, already mentioned,  is that

(3a)     (1/(2  π))(∫_(π/
              32)^(2  π - π/32)) ((φ(ξ))^2) dξ ≈ 0

This can be rewritten as

(3b)     h·h – Kh·h ≈ 0

where the matrix (Subscript[K, m, n]) (here denoted by K[[m,n]]) is given by:

K = 
        Table[If[m == n, N[(1/32)], 
            N[Sin[(m - n) π/32]/((m - n) π)]], {m, 1, 
            511}, {n, 1, 511}];

Indeed:

h . h - (K . h) . h

A final requirement is that

(4a) \(TraditionalForm\`\(\(\ \ \)\(\(\(\ \)\(∑\)\)\+\(\(\ \ \ \ \)\(k = 0\)\)\%\(\(\ \ \ \)\(63\)\)\)\)\)  ((φ(ξ+kπ/32))^2)≈ 1

Condition (3a) implies condition (1a), and in view of (0) it is almost equivalent to

(4b)    ((φ(ξ))^2) + ((φ(π/32 – ξ))^2)≈ 1   if    ξ ϵ [0, π/32].

Indeed:

Plot[φ[ξ]^2 + φ[π/32 -  ξ]^2 - 
          1, {ξ, 0, π/32}, 
        Ticks -> {{π/128, π/64, 3  π/128, π/32}, 
            Automatic}];

Condition (4a) is also equivalent to

(2b)    \(TraditionalForm\`\(\(\ \)\(∑\)\)\+\(n = 1\)\%511\) (Subscript[h, n]^(2 ))= h·h ≈ (1/64)

together with

(4c)    \(TraditionalForm\`\(\(\ \)\(∑\)\)\+\(n\  = \ 32  l\  + \ 1\)\%\(511  –\ 32\ l\)\)   (Subscript[h, n - 32  l]) (Subscript[h, n + 32  l ])≈ 0          (l =1, 2, ⋯, 7)

Indeed:

Table[{l, 
          Sum[h[[n - 32  l]] h[[n + 32  l]], {n, 1 + 32  l, 
              511 - 32  l}]}, {l, 1, 7}];
TableForm[%]

We still have to add a parity condition:

(5a)    φ(0) > 0

or equivalently

(5b)         \(TraditionalForm\`\(\(\ \)\(∑\)\)\+\(\(\ \ \)\(n = 1\)\)\%\(\(\ \ \)\(511\)\)\)  (Subscript[h, n ])> 0

The matrix   K   is symmetric and has eigenvalues as follows:

Chop[Eigenvalues[K]]

Hence
(3b)     h·h – Kh·h ≈ 0           together with               (1a)     (Subscript[h, n])= (Subscript[h, 512 - n])
will be satisfied if   h  is a linear combination of symmetric eigenvectors of  K   with eigenvalues close to 1.

There is a Toeplitz matrix  T   which commutes with   K :

T = 
        Table[Which[j == i - 1, j (511 - j)/2, j == i + 1, 
            i (511 - i)/2, j == i, 
            N[(((256 - j))^2) Cos[π/32]], True, 0], {i, 1, 511}, {j, 
            1, 511}];

KT = K . T - T . K;
Sum[Abs[KT[[i, j]]], {i, 1, 511}, {j, 1, 511}]

We pass to 256 × 256 versions K0   and   T0   of   K   resp.   T   (still symmetric matrices) which act on vectors with 256 coordinates which are obtained from symmetric vectors with 511 coordinates.

K0 = 
        Table[Which[
            1 ≤ m ≤ 255 && 
              1 ≤ n ≤ 255, 
            K[[m, n]] + K[[m, 512 - n]], 
            m == 256 && 1 ≤ n ≤ 255, 
            K[[256, n]] Sqrt[2], 
            1 ≤ m ≤ 255 && n == 256, 
            K[[m, 256]] Sqrt[2], True, K[[256, 256]]], {m, 1, 
            256}, {n, 1, 256}];

T0 = 
        Table[Which[
            1 ≤ m ≤ 255 && 
              1 ≤ n ≤ 255, 
            T[[m, n]] + T[[m, 512 - n]], 
            m == 256 && 1 ≤ n ≤ 255, 
            T[[256, n]] Sqrt[2], 
            1 ≤ m ≤ 255 && n == 256, 
            T[[m, 256]] Sqrt[2], True, T[[256, 256]]], {m, 1, 
            256}, {n, 1, 256}];

K0T0 = K0 . T0 - T0 . K0;

Sum[Abs[K0T0[[i, j]]], {i, 1, 256}, {j, 1, 256}]

Chop[Eigenvalues[K0]]

Table[1 - Eigenvalues[K0][[m]], {m, 1, 5}]

The five eigenvectors of   T0    corresponding to its highest eigenvalues are also eigenvectors of   K0    corresponding to its highest eigenvalues:

v0 = Table[Eigenvectors[T0][[m]], {m, 1, 5}];

Table[{1 - (K0 . v0[[m]])[[256]]/v0[[m, 256]], 
          Norm[K0 . 
                v0[[m]] - (K0 . v0[[m]])[[256]]/v0[[m, 
                        256]] v0[[m]]]}, {m, 1, 
          5}];
TableForm[%]

We extend these five eigenvectors   v0[[m]]   to symmetric eigenvectors   v[[m]]   of    K   and   T :

v = 
        Table[Which[n < 256, v0[[i, n]]/Sqrt[2], n > 256, 
            v0[[i, 512 - n]]/Sqrt[2], True, v0[[i, 256]]], {i, 1, 
            5}, {n, 1, 511}];

The five vectors   v[[m]]  form an orthonormal system:

Table[
        Chop[v[[m]] . v[[n]]], {m, 1, 5}, {n, 1, 
          5}];
TableForm[%]

Now we look for coefficients  a[i]  such that
h1 =\(TraditionalForm\`\(\(\ \)\(∑\)\)\+\(\(\ \ \)\(i = 1\)\)\%\(\(\ \ \)\(5\)\)\)  a[ i ] v[[ i ]]
satisfies     h1·h1 – Kh1·h1 ≈ 0   in the sense that
\(TraditionalForm\`\(\(\ \)\(∑\)\)\+\(\(\ \ \)\(l = 0\)\)\%\(\(\ \ \)\(7\)\)\) (((\(\(\(\ \)\(∑\)\)\+\(\(\ \ \)\(m = 
                          32  l + 1\)\)\%\(\(\ \ \)\(511  –\ 32  
                          l\)\)\)  h1[[ m – 32  l ]]  h1[[ m  +  32  l ]] – (1/64)(Subscript[δ, l, 0])))^2)     is minimal.

With the aid of the quadratic form

q = 
        Table[Sum[
            v[[i, m - 32  l]] v[[j, m + 32  l]], {m, 1 + 32  l, 
              511 - 32  l}], {l, 1, 7}, {i, 1, 5}, {j, 1, 5}];

we have to dind the following minimum:

FindMinimum[
      Sum[Sum[q[[n, i, j]] a[i] a[j], {i, 1, 5}, {j, 1, 5}]^2, {n, 1, 
            7}] + (Sum[a[i]^2, {i, 1, 5}] - 1/64)^2, {{a[1], 0.1}, {a[2], 
          0}, {a[3], 0}, {a[4], 0}, {a[5], 0}}]

Then h1 below is the desired vector:

h1 = 
        Sum[a[j] v[[j]], {j, 1, 
              5}] /. %[[2]];
If[Sum[h1[[n]], {n, 1, 511}] < 0, h1 = -h1];

We compare h1 with h

pl1 = 
        ListPlot[h1, PlotRange -> {All, All}, 
          DisplayFunction -> Identity];
pl = 
        ListPlot[h, PlotRange -> {All, All}, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{pl1, pl}], ImageSize -> 600];

Table[{l, 
          Sum[h1[[n - 32  l]] h1[[n + 32  l]], {n, 1 + 32  l, 
              511 - 32  l}], 
          Sum[h[[n - 32  l]] h[[n + 32  l]], {n, 1 + 32  l, 
              511 - 32  l}]}, {l, 1, 7}];
TableForm[%]

{h1 . h1 - 1/64, h . h - 1/64};
TableForm[%]

φ1[ξ_] := 
      h1[[256]] + 
        2  Sum[h1[[256 - n]] Cos[n ξ], {n, 1, 255}]

pl1 = 
        Plot[{φ1[ξ], UnitStep[π/64 - ξ]}, {ξ, 
            0, π/32}, PlotStyle -> {Black, Red}, 
          Ticks -> {{π/128, π/64, 3  π/128, π/32}, 
              Automatic}, 
          DisplayFunction -> Identity];
pl = 
        Plot[{φ[ξ], UnitStep[π/64 - ξ]}, {ξ, 
            0, π/32}, PlotStyle -> {Black, Red}, 
          Ticks -> {{π/128, π/64, 3  π/128, π/32}, 
              Automatic}, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{pl1, pl}], ImageSize -> 550];

pl1 = 
        Plot[φ1[ξ], {ξ, π/32, π}, 
          Ticks -> {{π/32, π}, Automatic}, 
          DisplayFunction -> Identity];
pl = 
        Plot[φ[ξ], {ξ, π/32, π}, 
          DisplayFunction -> Identity, 
          Ticks -> {{π/32, π}, Automatic}, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{pl1, pl}], ImageSize -> 550];

pl1 = 
        Plot[φ1[ξ]^2 + φ1[π/32 -  ξ]^2 - 
            1, {ξ, 0, π/32}, 
          Ticks -> {{π/128, π/64, 3  π/128, π/32}, 
              Automatic}, 
          DisplayFunction -> Identity];
pl = 
        Plot[φ[ξ]^2 + φ[π/32 -  ξ]^2 - 
            1, {ξ, 0, π/32}, 
          Ticks -> {{π/128, π/64, 3  π/128, π/32}, 
              Automatic}, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{pl1, pl}], ImageSize -> 550];

pl1 = 
        Plot[1 - φ1[ξ], {ξ, 0, π/256}, 
          Ticks -> {{π/512, π/256}, Automatic}, 
          DisplayFunction -> Identity];
pl = 
        Plot[1 - φ[ξ], {ξ, 0, π/256}, 
          Ticks -> {{π/512, π/256}, Automatic}, 
          DisplayFunction -> Identity];
Show[GraphicsArray[{pl1, pl}], ImageSize -> 550];

Norm[h1 - h]

ListPlot[h1 - h];