ref: 4e4db8fd09db5fc54df1572fab1eb93232d75ecd
parent: 7218a569ea0124a0a25dffd3e93b91bde67bdca7
author: robs <robs>
date: Sun Mar 18 11:32:50 EDT 2007
Changed deemph to be a true biquad for better accuracy.
--- a/ChangeLog
+++ b/ChangeLog
@@ -19,6 +19,7 @@
o Added soft-knee companding. (robs)
o Show (with --octave) compand transfer function. (robs)
o Allow e.g. "vol 6dB" (as well as "vol 6 dB"). (robs)
+ o Changed deemph to be a true biquad for better accuracy. (robs)
Other new features:
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -22,7 +22,7 @@
effects = band.h biquad.c biquad.h biquads.c chorus.c compand.c \
compandt.c compandt.h dcshift.c\
- deemph.h dither.c earwax.c echo.c echos.c fade.c FFT.c FFT.h filter.c\
+ dither.c earwax.c echo.c echos.c fade.c FFT.c FFT.h filter.c\
flanger.c mcompand.c mixer.c noiseprof.c \
noisered.c noisered.h pad.c pan.c phaser.c pitch.c polyphas.c \
rabbit.c rate.c repeat.c resample.c reverb.c reverse.c silence.c \
--- a/src/biquads.c
+++ b/src/biquads.c
@@ -130,7 +130,11 @@
static int deemph_getopts(eff_t effp, int n, char **argv) {- return sox_biquad_getopts(effp, n, argv, 0, 0, 0, 1, 2, "", filter_deemph);
+ biquad_t p = (biquad_t) effp->priv;
+ p->fc = 5283;
+ p->width = 0.4845;
+ p->gain = -9.477;
+ return sox_biquad_getopts(effp, n, argv, 0, 0, 0, 1, 2, "s", filter_deemph);
}
@@ -248,6 +252,13 @@
p->a2 = (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha;
break;
+ case filter_deemph: /* See deemph.plt for documentation */
+ if (effp->ininfo.rate != 44100) {+ sox_fail("Sample rate must be 44100 (audio-CD)");+ return SOX_EOF;
+ }
+ /* Falls through... */
+
case filter_highShelf: /* H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1)/(s^2 + (sqrt(A)/Q)*s + A) */
if (!A)
return SOX_EFF_NULL;
@@ -295,11 +306,6 @@
p->a1 = -2 * cos(w0);
p->a2 = 1 - sin(w0);
break;
-
- case filter_deemph: {- #include "deemph.h" /* Has own documentation */
- break;
- }
}
return sox_biquad_start(effp);
}
--- a/src/deemph.h
+++ /dev/null
@@ -1,109 +1,0 @@
-/*
- * July 5, 1991
- *
- * Deemphases Filter
- *
- * Fixed deemphasis filter for processing pre-emphasized audio cd samples
- * 09/02/98 (c) Heiko Eissfeldt
- * License: LGPL (Lesser Gnu Public License)
- *
- * This implements the inverse filter of the optional pre-emphasis stage as
- * defined by ISO 908 (describing the audio cd format).
- *
- * Background:
- * In the early days of audio cds, there were recording problems
- * with noise (for example in classical recordings). The high dynamics
- * of audio cds exposed these recording errors a lot.
- *
- * The commonly used solution at that time was to 'pre-emphasize' the
- * trebles to have a better signal-noise-ratio. That is trebles were
- * amplified before recording, so that they would give a stronger
- * signal compared to the underlying (tape)noise.
- *
- * For that purpose the audio signal was prefiltered with the following
- * frequency response (simple first order filter):
- *
- * V (in dB)
- * ^
- * |
- * | _________________
- * | /
- * | / |
- * | 20 dB / decade ->/ |
- * | / |
- * |____________________/_ _ |_ _ _ _ _ _ _ _ _ _ _ _ _ lg f
- * |0 dB | |
- * | | |
- * | | |
- * 3.1KHz ca. 10KHz
- *
- * So the recorded audio signal has amplified trebles compared to the
- * original.
- * HiFi cd players do correct this by applying an inverse filter
- * automatically, the cd-rom drives or cd burners used by digital
- * sampling programs (like cdda2wav) however do not.
- *
- * So, this is what this effect does.
- *
- * Here is the gnuplot file for the frequency response
- of the deemphasis. The error is below +-0.1dB
-
--------- Start of gnuplot file ---------------------
-# first define the ideal filter. We use the tenfold sampling frequency.
-T=1./441000.
-OmegaU=1./15E-6
-OmegaL=15./50.*OmegaU
-V0=OmegaL/OmegaU
-H0=V0-1.
-B=V0*tan(OmegaU*T/2.)
-# the coefficients follow
-a1=(B - 1.)/(B + 1.)
-b0=(1.0 + (1.0 - a1) * H0/2.)
-b1=(a1 + (a1 - 1.0) * H0/2.)
-# helper variables
-D=b1/b0
-o=2*pi*T
-H2(f)=b0*sqrt((1+2*cos(f*o)*D+D*D)/(1+2*cos(f*o)*a1+a1*a1))
-#
-# now approximate the ideal curve with a fitted one for sampling
-frequency
-# of 44100 Hz. Fitting parameters are
-# amplification at high frequencies V02
-# and tau of the upper edge frequency OmegaU2 = 2 *pi * f(upper)
-T2=1./44100.
-V02=0.3365
-OmegaU2=1./19E-6
-B2=V02*tan(OmegaU2*T2/2.)
-# the coefficients follow
-a12=(B2 - 1.)/(B2 + 1.)
-b02=(1.0 + (1.0 - a12) * (V02-1.)/2.)
-b12=(a12 + (a12 - 1.0) * (V02-1.)/2.)
-# helper variables
-D2=b12/b02
-o2=2*pi*T2
-H(f)=b02*sqrt((1+2*cos(f*o2)*D2+D2*D2)/(1+2*cos(f*o2)*a12+a12*a12))
-# plot best, real, ideal, level with halved attenuation,
-# level at full attentuation, 10fold magnified error
-set logscale x
-set grid xtics ytics mxtics mytics
-plot [f=1000:20000] [-12:2] 20*log10(H(f)),20*log10(H2(f)),
-20*log10(OmegaL/(2*
-pi*f)), 0.5*20*log10(V0), 20*log10(V0), 200*log10(H(f)/H2(f))
-pause -1 "Hit return to continue"
--------- End of gnuplot file ---------------------
-*/
-
-
-/* filter coefficients */
-p->a1 = -0.62786881719628784282;
-p->b0 = 0.45995451989513153057;
-p->b1 = -0.08782333709141937339;
-
-
-/* The sample-rate must be 44100 as this has been hard-coded into the
- * pre-calculated filter coefficients.
- */
-if (effp->ininfo.rate != 44100) {- sox_fail("Sample rate must be 44100 (audio-CD)");- return SOX_EOF;
-}
--- /dev/null
+++ b/src/deemph.plt
@@ -1,0 +1,121 @@
+# 15/50us EIAJ de-emphasis filter for CD/DAT
+#
+# 09/02/98 (c) Heiko Eissfeldt
+#
+# 18/03/07 robs@users.sourceforge.net: changed to biquad for better accuracy.
+#
+# License: LGPL (Lesser Gnu Public License)
+#
+# This implements the inverse filter of the optional pre-emphasis stage
+# as defined by ISO 908 (describing the audio cd format).
+#
+# Background: In the early days of audio cds, there were recording
+# problems with noise (for example in classical recordings). The high
+# dynamics of audio cds exposed these recording errors a lot.
+#
+# The commonly used solution at that time was to 'pre-emphasize' the
+# trebles to have a better signal-noise-ratio. That is trebles were
+# amplified before recording, so that they would give a stronger signal
+# compared to the underlying (tape) noise.
+#
+# For that purpose the audio signal was prefiltered with the following
+# frequency response (simple first order filter):
+#
+# V (in dB)
+# ^
+# |
+# |~10dB _________________
+# | /
+# | / |
+# | 20dB / decade ->/ |
+# | / |
+# |____________________/_ _ |_ _ _ _ _ _ _ _ _ Frequency
+# |0 dB | |
+# | | |
+# | | |
+# 3.1kHz ~10kHz
+#
+# So the recorded audio signal has amplified trebles compared to the
+# original. HiFi cd players do correct this by applying an inverse
+# filter automatically, the cd-rom drives or cd burners used by digital
+# sampling programs (like cdda2wav) however do not.
+#
+# So, this is what this effect does.
+#
+# This is the gnuplot file for the frequency response of the deemphasis.
+#
+# The absolute error is <=0.04dB up to ~12kHz, and <=0.06dB up to 20kHz.
+
+# First define the ideal filter:
+
+# Filter parameters
+T=1./441000. # we use the tenfold sampling frequency
+OmegaU=1./15e-6
+OmegaL=15./50.*OmegaU
+
+# Calculate filter coefficients
+V0=OmegaL/OmegaU
+H0=V0 - 1.
+B=V0*tan(OmegaU*T/2.)
+a1=(B - 1.)/(B + 1.)
+b0=(1.0 + (1.0 - a1) * H0/2.)
+b1=(a1 + (a1 - 1.0) * H0/2.)
+
+# helper variables
+D=b1/b0
+o=2*pi*T
+
+# Ideal transfer function
+Hi(f)=b0*sqrt((1 + 2*cos(f*o)*D + D*D)/(1 + 2*cos(f*o)*a1 + a1*a1))
+
+# Now use a biquad (RBJ high shelf) with sampling frequency of 44100 Hz
+# to approximate the ideal curve:
+
+# Filter parameters
+m_t=1./44100.
+m_gain=-9.477
+m_slope=.4845
+m_f0=5283
+
+# Calculate filter coefficients
+m_A=exp(m_gain/40.*log(10.))
+m_w0=2.*pi*m_f0*m_t
+m_alpha=sin(m_w0)/2.*sqrt((m_A+1./m_A)*(1./m_slope-1.)+2.)
+m_b0=m_A*((m_A+1.)+(m_A-1.)*cos(m_w0)+2.*sqrt(m_A)*m_alpha)
+m_b1=-2.*m_A*((m_A-1.)+(m_A+1.)*cos(m_w0))
+m_b2=m_A*((m_A+1.)+(m_A-1.)*cos(m_w0)-2.*sqrt(m_A)*m_alpha)
+m_a0=(m_A+1.)-(m_A-1.)*cos(m_w0)+2.*sqrt(m_A)*m_alpha
+m_a1=2.*((m_A-1.)-(m_A+1.)*cos(m_w0))
+m_a2=(m_A+1.)-(m_A-1.)*cos(m_w0)-2.*sqrt(m_A)*m_alpha
+m_b2=m_b2/m_a0
+m_b1=m_b1/m_a0
+m_b0=m_b0/m_a0
+m_a2=m_a2/m_a0
+m_a1=m_a1/m_a0
+
+# helper variables
+m_o=2*pi*m_t
+
+# Best fit transfer function
+Hb(f)=sqrt(\
+ (m_b0*m_b0 + m_b1*m_b1 + m_b2*m_b2 +\
+ 2.*(m_b0*m_b1 + m_b1*m_b2)*cos(f*m_o) +\
+ 2.*(m_b0*m_b2)* cos(2.*f*m_o)) /\
+ (1. + m_a1*m_a1 + m_a2*m_a2 +\
+ 2.*(m_a1 + m_a1*m_a2)*cos(f*m_o) +\
+ 2.*m_a2*cos(2.*f*m_o)))
+
+# plot real, best, ideal, level with halved attenuation,
+# level at full attentuation, 10fold magnified error
+set logscale x
+set grid xtics ytics mxtics mytics
+set key left bottom
+plot [f=1000:20000] [-12:2] \
+20*log10(Hi(f)),\
+20*log10(Hb(f)),\
+20*log10(OmegaL/(2* pi*f)),\
+0.5*20*log10(V0),\
+20*log10(V0),\
+200*log10(Hb(f)/Hi(f))
+
+pause -1 "Hit return to continue"