1# 15/50us EIAJ de-emphasis filter for CD/DAT
2#
3# 09/02/98 (c) Heiko Eissfeldt
4#
5# 18/03/07 robs@users.sourceforge.net: changed to biquad for slightly
6# better accuracy.
7#
8# License: LGPL (Lesser Gnu Public License)
9#
10# This implements the inverse filter of the optional pre-emphasis stage
11# as defined by IEC 60908 (describing the audio cd format).
12#
13# Background: In the early days of audio cds, there were recording
14# problems with noise (for example in classical recordings). The high
15# dynamics of audio cds exposed these recording errors a lot.
16#
17# The commonly used solution at that time was to 'pre-emphasize' the
18# trebles to have a better signal-noise-ratio. That is trebles were
19# amplified before recording, so that they would give a stronger signal
20# compared to the underlying (tape) noise.
21#
22# For that purpose the audio signal was prefiltered with the following
23# frequency response (simple first order filter):
24#
25# V (in dB)
26# ^
27# |
28# |~10dB                    _________________
29# |                        /
30# |                       / |
31# |      20dB / decade ->/  |
32# |                     /   |
33# |____________________/_ _ |_ _ _ _ _ _ _ _ _ Frequency
34# |0 dB                |    |
35# |                    |    |
36# |                    |    |
37#                 3.1kHz    ~10kHz
38#
39# So the recorded audio signal has amplified trebles compared to the
40# original.  HiFi cd players do correct this by applying an inverse
41# filter automatically, the cd-rom drives or cd burners used by digital
42# sampling programs (like cdda2wav) however do not.
43#
44# So, this is what this effect does.
45#
46# This is the gnuplot file for the frequency response of the deemphasis.
47#
48# The absolute error is <=0.04dB up to ~12kHz, and <=0.06dB up to 20kHz.
49
50# First define the ideal filter:
51
52# Filter parameters
53T = 1. / 441000.          # we use the tenfold sampling frequency
54OmegaU = 1. / 15e-6
55OmegaL = 15. / 50. * OmegaU
56
57# Calculate filter coefficients
58V0 = OmegaL / OmegaU
59H0 = V0 - 1.
60B = V0 * tan(OmegaU * T / 2.)
61A1 = (B - 1.) / (B + 1.)
62B0 = (1. + (1. - A1) * H0 / 2.)
63B1 = (A1 + (A1 - 1.) * H0 / 2.)
64
65# helper variables
66D = B1 / B0
67O = 2 * pi * T
68
69# Ideal transfer function
70Hi(f) = B0*sqrt((1 + 2*cos(f*O)*D + D*D)/(1 + 2*cos(f*O)*A1 + A1*A1))
71
72# Now use a biquad (RBJ high shelf) with sampling frequency of 44100Hz
73# to approximate the ideal curve:
74
75# Filter parameters
76t = 1. / 44100.
77gain = -9.477
78slope = .4845
79f0 = 5283
80
81# Calculate filter coefficients
82A = exp(gain / 40. * log(10.))
83w0 = 2. * pi * f0 * t
84alpha = sin(w0) / 2. * sqrt((A + 1. / A) * (1. / slope - 1.) + 2.)
85b0 = A * ((A + 1.) + (A - 1.) * cos(w0) + 2. * sqrt(A) * alpha)
86b1 = -2. * A * ((A - 1.) + (A + 1.) * cos(w0))
87b2 = A * ((A + 1.) + (A - 1.) * cos(w0) - 2. * sqrt(A) * alpha)
88a0 = (A + 1.) - (A - 1.) * cos(w0) + 2. * sqrt(A) * alpha
89a1 = 2. * ((A - 1.) - (A + 1.) * cos(w0))
90a2 = (A + 1.) - (A - 1.) * cos(w0) - 2. * sqrt(A) * alpha
91b2 = b2 / a0
92b1 = b1 / a0
93b0 = b0 / a0
94a2 = a2 / a0
95a1 = a1 / a0
96
97# helper variables
98o = 2 * pi * t
99
100# Best fit transfer function
101Hb(f) = sqrt((b0*b0 + b1*b1 + b2*b2 +\
102    2.*(b0*b1 + b1*b2)*cos(f*o) + 2.*(b0*b2)* cos(2.*f*o)) /\
103  (1. + a1*a1 + a2*a2 + 2.*(a1 + a1*a2)*cos(f*o) + 2.*a2*cos(2.*f*o)))
104
105# plot real, best, ideal, level with halved attenuation,
106#   level at full attentuation, 10fold magnified error
107set logscale x
108set grid xtics ytics mxtics mytics
109set key left bottom
110plot [f=1000:20000] [-12:2] \
11120 * log10(Hi(f)),\
11220 * log10(Hb(f)),\
11320 * log10(OmegaL/(2 * pi * f)),\
114.5 * 20 * log10(V0),\
11520 * log10(V0),\
116200 * log10(Hb(f)/Hi(f))
117
118pause -1 "Hit return to continue"
119