Use GNU Radio emphasis filters, add filters generation scripts
This commit is contained in:
parent
eef364b1a9
commit
98cd2404f5
|
@ -25,7 +25,7 @@
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
const float DEEMPHASIS_GAIN_DB = 0.0F;
|
const float DEEMPHASIS_GAIN_DB = 0.0F;
|
||||||
const float PREEMPHASIS_GAIN_DB = 30.0F;
|
const float PREEMPHASIS_GAIN_DB = 0.0F;
|
||||||
const float FILTER_GAIN_DB = 0.0F;
|
const float FILTER_GAIN_DB = 0.0F;
|
||||||
const unsigned int FM_MASK = 0x00000FFFU;
|
const unsigned int FM_MASK = 0x00000FFFU;
|
||||||
|
|
||||||
|
@ -39,8 +39,8 @@ m_filterStage1(NULL),
|
||||||
m_filterStage2(NULL),
|
m_filterStage2(NULL),
|
||||||
m_filterStage3(NULL)
|
m_filterStage3(NULL)
|
||||||
{
|
{
|
||||||
m_preemphasis = new CIIRDirectForm1Filter(0.38897032f, -0.32900053f, 0.0f, 1.0f, 0.28202918f, 0.0f, PREEMPHASIS_GAIN_DB);
|
m_preemphasis = new CIIRDirectForm1Filter(8.315375384336983F,-7.03334621603483F,0.0F,1.0F,0.282029168302153F,0.0F, PREEMPHASIS_GAIN_DB);
|
||||||
m_deemphasis = new CIIRDirectForm1Filter(1.0f,0.28202918f, 0.0f, 0.38897032f, -0.32900053f, 0.0f, DEEMPHASIS_GAIN_DB);
|
m_deemphasis = new CIIRDirectForm1Filter(0.07708787090460224F,0.07708787090460224F,0.0F,1.0F,-0.8458242581907955F,0.0F, DEEMPHASIS_GAIN_DB);
|
||||||
|
|
||||||
//cheby type 1 0.2dB cheby type 1 3rd order 300-2700Hz fs=8000
|
//cheby type 1 0.2dB cheby type 1 3rd order 300-2700Hz fs=8000
|
||||||
m_filterStage1 = new CIIRDirectForm1Filter(0.29495028f, 0.0f, -0.29495028f, 1.0f, -0.61384624f, -0.057158668f, FILTER_GAIN_DB);
|
m_filterStage1 = new CIIRDirectForm1Filter(0.29495028f, 0.0f, -0.29495028f, 1.0f, -0.61384624f, -0.057158668f, FILTER_GAIN_DB);
|
||||||
|
|
43
Tools/DeEmphasis.py
Normal file
43
Tools/DeEmphasis.py
Normal file
|
@ -0,0 +1,43 @@
|
||||||
|
#based on https://github.com/gnuradio/gnuradio/blob/master/gr-analog/python/analog/fm_emph.py
|
||||||
|
|
||||||
|
import math
|
||||||
|
import cmath
|
||||||
|
import numpy as np
|
||||||
|
import scipy.signal as signal
|
||||||
|
import pylab as pl
|
||||||
|
|
||||||
|
tau = 750e-6
|
||||||
|
fs = 8000
|
||||||
|
fh = 2700
|
||||||
|
|
||||||
|
# Digital corner frequency
|
||||||
|
w_c = 1.0 / tau
|
||||||
|
|
||||||
|
# Prewarped analog corner frequency
|
||||||
|
w_ca = 2.0 * fs * math.tan(w_c / (2.0 * fs))
|
||||||
|
|
||||||
|
# Resulting digital pole, zero, and gain term from the bilinear
|
||||||
|
# transformation of H(s) = w_ca / (s + w_ca) to
|
||||||
|
# H(z) = b0 (1 - z1 z^-1)/(1 - p1 z^-1)
|
||||||
|
k = -w_ca / (2.0 * fs)
|
||||||
|
z1 = -1.0
|
||||||
|
p1 = (1.0 + k) / (1.0 - k)
|
||||||
|
b0 = -k / (1.0 - k)
|
||||||
|
|
||||||
|
btaps = [ b0 * 1.0, b0 * -z1, 0 ]
|
||||||
|
ataps = [ 1.0, -p1, 0 ]
|
||||||
|
|
||||||
|
# Since H(s = 0) = 1.0, then H(z = 1) = 1.0 and has 0 dB gain at DC
|
||||||
|
|
||||||
|
|
||||||
|
taps = np.concatenate((btaps, ataps), axis=0)
|
||||||
|
print("Taps")
|
||||||
|
print(*taps, "", sep=",", end="\n")
|
||||||
|
|
||||||
|
f,h = signal.freqz(btaps,ataps, fs=fs)
|
||||||
|
pl.plot(f, 20*np.log10(np.abs(h)))
|
||||||
|
pl.xlabel('frequency/Hz')
|
||||||
|
pl.ylabel('gain/dB')
|
||||||
|
pl.ylim(top=0,bottom=-30)
|
||||||
|
pl.xlim(left=0, right=fh*2.5)
|
||||||
|
pl.show()
|
51
Tools/PreEmphasis.py
Normal file
51
Tools/PreEmphasis.py
Normal file
|
@ -0,0 +1,51 @@
|
||||||
|
#based on https://github.com/gnuradio/gnuradio/blob/master/gr-analog/python/analog/fm_emph.py
|
||||||
|
|
||||||
|
import math
|
||||||
|
import cmath
|
||||||
|
import numpy as np
|
||||||
|
import scipy.signal as signal
|
||||||
|
import pylab as pl
|
||||||
|
|
||||||
|
tau = 750e-6
|
||||||
|
fs = 8000
|
||||||
|
fh = 2700
|
||||||
|
|
||||||
|
# Digital corner frequencies
|
||||||
|
w_cl = 1.0 / tau
|
||||||
|
w_ch = 2.0 * math.pi * fh
|
||||||
|
|
||||||
|
# Prewarped analog corner frequencies
|
||||||
|
w_cla = 2.0 * fs * math.tan(w_cl / (2.0 * fs))
|
||||||
|
w_cha = 2.0 * fs * math.tan(w_ch / (2.0 * fs))
|
||||||
|
|
||||||
|
# Resulting digital pole, zero, and gain term from the bilinear
|
||||||
|
# transformation of H(s) = (s + w_cla) / (s + w_cha) to
|
||||||
|
# H(z) = b0 (1 - z1 z^-1)/(1 - p1 z^-1)
|
||||||
|
kl = -w_cla / (2.0 * fs)
|
||||||
|
kh = -w_cha / (2.0 * fs)
|
||||||
|
z1 = (1.0 + kl) / (1.0 - kl)
|
||||||
|
p1 = (1.0 + kh) / (1.0 - kh)
|
||||||
|
b0 = (1.0 - kl) / (1.0 - kh)
|
||||||
|
|
||||||
|
# Since H(s = infinity) = 1.0, then H(z = -1) = 1.0 and
|
||||||
|
# this filter has 0 dB gain at fs/2.0.
|
||||||
|
# That isn't what users are going to expect, so adjust with a
|
||||||
|
# gain, g, so that H(z = 1) = 1.0 for 0 dB gain at DC.
|
||||||
|
w_0dB = 2.0 * math.pi * 0.0
|
||||||
|
g = abs(1.0 - p1 * cmath.rect(1.0, -w_0dB)) \
|
||||||
|
/ (b0 * abs(1.0 - z1 * cmath.rect(1.0, -w_0dB)))
|
||||||
|
|
||||||
|
btaps = [ g * b0 * 1.0, g * b0 * -z1, 0]
|
||||||
|
ataps = [ 1.0, -p1, 0]
|
||||||
|
|
||||||
|
taps = np.concatenate((btaps, ataps), axis=0)
|
||||||
|
print("Taps")
|
||||||
|
print(*taps, "", sep=",", end="\n")
|
||||||
|
|
||||||
|
f,h = signal.freqz(btaps,ataps, fs=fs)
|
||||||
|
pl.plot(f, 20*np.log10(np.abs(h)))
|
||||||
|
pl.xlabel('frequency/Hz')
|
||||||
|
pl.ylabel('gain/dB')
|
||||||
|
pl.ylim(top=30,bottom=0)
|
||||||
|
pl.xlim(left=0, right=fh*2.5)
|
||||||
|
pl.show()
|
Loading…
Reference in a new issue