Designing De-emphasis Digital Filter for Old CDs Part 2 : Frequency Sampling Method

On the previous blog post, we've got a 6th degree polynomial function of de-emphasis curve Equation 1:

$f\mathrm{dB}(x)\; =\; -0.0000029212346025337816x^6+0.00020291497408909238x^5-0.0054888099286801205x^4+0.071110615465301924x^3-0.40078216359333169x^2-0.11354738870338571x$

$x\; :\; frequency\; (kHz)$
$f_{\mathrm{dB}}(x):\; filter\; gain\; (dB)$

Note: Math equations of this article uses MathML and it seems it is displayed correctly only on Firefox (as of November 2020).

This decibel is root power quantity therefore actual gain value f(x) is:

Next step is to create filter gain table of equal spacing frequency using this function.

And create the filter, evaluate its frequency response error from original de-emphasis table and choose the best FIR filter tap number M.

 

Example : Calculation of M=9 taps FIR filter coeffs

When sampling frequency==44100 Hz and Desired FIR filter taps M==9,

Frequency sampling index $k=0,1,2,\; ..,\frac{M-1}{2}=4$

$k=0,1,2,3\; and\; 4$

Frequency sampling angle frequency ${\omega }_{\mathrm{k}}=\frac{2\pi k}{M}:$

• $k=0: {\omega }_0=\; 0, {\mathrm{freq}}_0=\; 0\; Hz$
  
• $k=1: {\omega }_1=\; 2\pi /9, {\mathrm{freq}}_1=\; (44100/2\pi )(2\pi /9)\; =\; 44100\; /\; 9\; Hz\; =\; 4900\; Hz$
• $k=2: {\omega }_2=\; 4\pi /9, {\mathrm{freq}}_2=\; (44100/2\pi )(4\pi /9)\; =\; 44100\; *\; 2\; /\; 9\; Hz\; =\; 9800\; Hz$
• $k=3: {\omega }_3=\; 6\pi /9, {\mathrm{freq}}_3=\; (44100/2\pi )(6\pi /9)\; =\; 44100\; *\; 3\; /\; 9\; Hz\; =\; 14700\; Hz$
• $k=4: {\omega }_4=\; 8\pi /9, {\mathrm{freq}}_4=\; (44100/2\pi )(8\pi /9)\; =\; 44100\; *\; 4\; /\; 9\; Hz\; =\; 19600\; Hz$

Get filter gain on those frequencies $Hr({\omega }_{\mathrm{k}})$ using Equation (2):

• $Hr({\omega }_0)\; =\; f(0)\; =\; 1$
• $Hr({\omega }_1)\; =\; f(4.9)\; =\; 0.599479869$
• $Hr({\omega }_2)\; =\; f(9.8)\; =\; 0.419371436$
• $Hr({\omega }_3)\; =\; f(14.7)\; =\; 0.359695479$
• $Hr({\omega }_4)\; =\; f(19.6)\; =\; 0.33620803$

Then G(k) is calculated from $Hr({\omega }_{\mathrm{k}})$ by $G(k)\; ={\mathrm{(-1)}}^{\mathrm{k}}Hr({\omega }_{\mathrm{k}})\; \dots \dots (3)\; :$

• $G(0)\; ={\mathrm{(-1)}}^{0}Hr({\omega }_0)\; =\; 1$
• $G(1)\; ={\mathrm{(-1)}}^{1}Hr({\omega }_1)\; =\; -0.599479869$
• $G(2)\; ={\mathrm{(-1)}}^{2}Hr({\omega }_2)\; =\; 0.419371436$
• $G(3)\; ={\mathrm{(-1)}}^{3}Hr({\omega }_3)\; =\; -0.359695479$
• $G(4)\; ={\mathrm{(-1)}}^{4}Hr({\omega }_4)\; =\; 0.33620803$

FIR filter coefficients h(n) can be calculated by $h(n)\; =\frac{1}{M}\left\{G(0)+2\sum _{\mathrm{k=1}}^{\mathrm{(M-1)/2}}G(k)cos\frac{2\pi k(n+1/2)}{M}\right\}\dots (4)$

• $h(0)\; =\; 0.0303254491484693$
• $h(1)\; =\; 0.0404812914444444$
• $h(2)\; =\; 0.0639379770301665$
• $h(3)\; =\; 0.119171414154698$
• $h(4)\; =\; 0.492167736444444$

And finally, This FIR filter is symmetry shape, therefore h(8) = h(0), h(7) = h(1), h(6) = h(2), h(5) = h(3).

• $h(0)\; =\; 0.0303254491484693$
• $h(1)\; =\; 0.0404812914444444$
• $h(2)\; =\; 0.0639379770301665$
• $h(3)\; =\; 0.119171414154698$
• $h(4)\; =\; 0.492167736444444$
• $h(5)\; =\; 0.119171414154698$
• $h(6)\; =\; 0.0639379770301665$
• $h(7)\; =\; 0.0404812914444444$
• $h(8)\; =\; 0.0303254491484693$

We've got all the 9 FIR coefficient values h(n).

Evaluating Frequency Response of FIR filter

Now one FIR filter coefficients h(n) is available for testing. FIR filter gain of arbitrary angular frequency ω can be calculated using the following equation:

$Gain(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)e^{\mathrm{-jk\omega }}\dots \dots (5)$

Real part and imaginary part of (5) can be calculated separately:

${\mathrm{Gain}}_{\mathrm{real}}(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)cos(-k\omega )\; \dots \dots (5r)$
${\mathrm{Gain}}_{\mathrm{imaginary}}(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)sin(-k\omega )\; \dots \dots (5i)$

And FIR filter Gain magnitude is calculated as follows:

Finally Gain in decibel is: ${\mathrm{Gain}}_{\mathrm{dB}}(\omega )\; =\; 20{\log }_{10}\left\{{\mathrm{Gain}}_{\mathrm{magnitude}}(\omega )\right\}\dots \dots (7)$

For our M=9 FIR filter, frequency response can be calculated using equation (7):

Frequency(kHz)	${\mathrm{Gain}}_{\mathrm{dB}}(\omega )$
0	0
1	-0.214929748
2	-0.847460501
3	-1.855702469
4	-3.153141013
5	-4.58835657
6	-5.939107409
7	-6.960958392
8	-7.506314657
9	-7.624120111
10	-7.521831664
11	-7.433783042
12	-7.52486218
13	-7.862458657
14	-8.419180318
15	-9.081142015
16	-9.670162018
17	-10.00514702
18	-9.997371674
19	-9.707880014
20	-9.30607583
21	-8.974681539
22	-8.840811148

Table 1

Comparing to the original de-emphasis table, max error of this FIR filter is 0.87 dB on 7kHz, It is poor and may be M=9 is too small.

Finding Optimal FIR filter taps M

I'd like to have < 0.1 dB error of FIR filter with minimum filter taps. Calculated error from the original table on several M using DesignFrequencySamplingFilter and compared their performance.

Filter taps M	Max error from the original table
9	0.871 dB
17	0.217 dB
19	0.140 dB
21	0.172 dB
23	0.104 dB
25	0.118 dB
27	0.0779 dB
29	0.0802 dB
31	0.0795 dB
33	0.0770 dB

Table 2

From the table 2, M=27 is the most desirable one.

On M=27,  maximum error from the original table is 0.0779 dB on 2kHz.

Resulted linear-phase FIR filter coefficients h(n) of M=27 for 44.1kHz PCM is as follows:

0.00087829953598830856,
0.00073354073461322569,
0.0013059505528472161,
0.00089158366884379073,
0.0022743712962963354,
0.0017509721612062167,
0.0046856769010995523,
0.0049418323357026065,
0.011621996337245248,
0.017825153275235591,
0.034805918374128435,
0.057946349576219941,
0.11652885832464696,
0.48761899385185181,
0.11652885832464696,
0.057946349576219941,
0.034805918374128435,
0.017825153275235591,
0.011621996337245248,
0.0049418323357026065,
0.0046856769010995523,
0.0017509721612062167,
0.0022743712962963354,
0.00089158366884379073,
0.0013059505528472161,
0.00073354073461322569,
0.00087829953598830856

Fig.1 M=27 de-emphasis FIR Filter frequency response

This linear-phase de-emphasis FIR filter is available on DSP menu of PlayPcmWin 5.0.79. Source code is https://sourceforge.net/p/playpcmwin/code/HEAD/tree/PlayPcmWin/WasapiIODLL/WWAudioFilterDeEmphasis.cpp

 

Continued to Part 3 (Improving accuracy): https://yamamoto2002.blogspot.com/2020/11/designing-de-emphasis-digital-filter.html

 

 Reference

 J.G. Proakis & D.G. Manolakis: Digital Signal Processing, 4th edition, 2007, Chapter 10, pp. 671-678