src/deemph.plt

   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
  53 T = 1. / 441000.          # we use the tenfold sampling frequency
  54 OmegaU = 1. / 15e-6
  55 OmegaL = 15. / 50. * OmegaU
  56
  57 # Calculate filter coefficients
  58 V0 = OmegaL / OmegaU
  59 H0 = V0 - 1.
  60 B = V0 * tan(OmegaU * T / 2.)
  61 A1 = (B - 1.) / (B + 1.)
  62 B0 = (1. + (1. - A1) * H0 / 2.)
  63 B1 = (A1 + (A1 - 1.) * H0 / 2.)
  64
  65 # helper variables
  66 D = B1 / B0
  67 O = 2 * pi * T
  68
  69 # Ideal transfer function
  70 Hi(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
  76 t = 1. / 44100.
  77 gain = -9.477
  78 slope = .4845
  79 f0 = 5283
  80
  81 # Calculate filter coefficients
  82 A = exp(gain / 40. * log(10.))
  83 w0 = 2. * pi * f0 * t
  84 alpha = sin(w0) / 2. * sqrt((A + 1. / A) * (1. / slope - 1.) + 2.)
  85 b0 = A * ((A + 1.) + (A - 1.) * cos(w0) + 2. * sqrt(A) * alpha)
  86 b1 = -2. * A * ((A - 1.) + (A + 1.) * cos(w0))
  87 b2 = A * ((A + 1.) + (A - 1.) * cos(w0) - 2. * sqrt(A) * alpha)
  88 a0 = (A + 1.) - (A - 1.) * cos(w0) + 2. * sqrt(A) * alpha
  89 a1 = 2. * ((A - 1.) - (A + 1.) * cos(w0))
  90 a2 = (A + 1.) - (A - 1.) * cos(w0) - 2. * sqrt(A) * alpha
  91 b2 = b2 / a0
  92 b1 = b1 / a0
  93 b0 = b0 / a0
  94 a2 = a2 / a0
  95 a1 = a1 / a0
  96
  97 # helper variables
  98 o = 2 * pi * t
  99
 100 # Best fit transfer function
 101 Hb(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
 107 set logscale x
 108 set grid xtics ytics mxtics mytics
 109 set key left bottom
 110 plot [f=1000:20000] [-12:2] \
 111 20 * log10(Hi(f)),\
 112 20 * log10(Hb(f)),\
 113 20 * log10(OmegaL/(2 * pi * f)),\
 114 .5 * 20 * log10(V0),\
 115 20 * log10(V0),\
 116 200 * log10(Hb(f)/Hi(f))
 117
 118 pause -1 "Hit return to continue"