Merge pull request #479 from atsampson/deemphasis

Simplify the video deemphasis filter generation

lddecode/core.py

@@ -141,7 +141,7 @@ def calclinelen(SP, mult, mhz):
     'audio_notchwidth': 350000,
     'audio_notchorder': 2,
 
-    'video_deemp': (120*.32, 320*.32),
+    'video_deemp': (120e-9, 320e-9),
 
     # This BPF is similar but not *quite* identical to what Pioneer did
     'video_bpf_low': 3400000, 
@@ -172,7 +172,7 @@ def calclinelen(SP, mult, mhz):
     'audio_notchwidth': 350000,
     'audio_notchorder': 2,
 
-    'video_deemp': (120*.32, 320*.32),
+    'video_deemp': (120e-9, 320e-9),
 
     'video_bpf_low': 3800000, 
     'video_bpf_high': 12500000,
@@ -199,7 +199,7 @@ def calclinelen(SP, mult, mhz):
     'audio_notchwidth': 200000,
     'audio_notchorder': 2,
 
-    'video_deemp': (100*.30, 400*.30),
+    'video_deemp': (100e-9, 400e-9),
 
     # XXX: guessing here!
     'video_bpf_low': 2700000, 
@@ -227,7 +227,7 @@ def calclinelen(SP, mult, mhz):
     'audio_notchwidth': 200000,
     'audio_notchorder': 2,
 
-    'video_deemp': (100*.30, 400*.30),
+    'video_deemp': (100e-9, 400e-9),
 
     # XXX: guessing here!
     'video_bpf_low': 3200000, 
@@ -452,15 +452,13 @@ def computevideofilters(self):
         video_hpf = sps.butter(DP['video_hpf_order'], DP['video_hpf_freq']/self.freq_hz_half, 'high')
         SF['Fvideo_hpf'] = filtfft(video_hpf, self.blocklen)
 
-        # The deemphasis filter.  This math is probably still quite wrong, but with the right values it works
-        deemp0, deemp1 = DP['video_deemp']
-        [tf_b, tf_a] = sps.zpk2tf([-deemp1*(10**-10)], [-deemp0*(10**-10)], deemp0 / deemp1)
-        SF['Fdeemp'] = filtfft(sps.bilinear(tf_b, tf_a, 1.0/self.freq_hz_half), self.blocklen)
+        # The deemphasis filter
+        deemp1, deemp2 = DP['video_deemp']
+        SF['Fdeemp'] = filtfft(emphasis_iir(deemp1, deemp2, self.freq_hz), self.blocklen)
 
         # The direct opposite of the above, used in test signal generation
-        [tf_b, tf_a] = sps.zpk2tf([-deemp0*(10**-10)], [-deemp1*(10**-10)], deemp1 / deemp0)
-        SF['Femp'] = filtfft(sps.bilinear(tf_b, tf_a, 1.0/self.freq_hz_half), self.blocklen)
-
+        SF['Femp'] = filtfft(emphasis_iir(deemp2, deemp1, self.freq_hz), self.blocklen)
+
         # Post processing:  lowpass filter + deemp
         SF['FVideo'] = SF['Fvideo_lpf'] * SF['Fdeemp'] 
         #SF['FVideo'] = SF['Fdeemp'] 

lddecode/utils.py

@@ -521,6 +521,18 @@ def __call__(self, infile, sample, readlen):
 # This can be complex multiplied with the raw RF to do a good chunk of real demoduation work
 #fft_hilbert = np.fft.fft(hilbert_filter, blocklen)
 
+def emphasis_iir(t1, t2, fs):
+    """Generate an IIR filter for 6dB/octave pre-emphasis (t1 > t2) or
+    de-emphasis (t1 < t2), given time constants for the two corners."""
+
+    # Convert time constants to frequencies, and pre-warp for bilinear transform
+    w1 = 2 * fs * np.tan((1 / t1) / (2 * fs))
+    w2 = 2 * fs * np.tan((1 / t2) / (2 * fs))
+
+    # Zero at t1, pole at t2
+    tf_b, tf_a = sps.zpk2tf([-w1], [-w2], w2 / w1)
+    return sps.bilinear(tf_b, tf_a, fs)
+
 # This converts a regular B, A filter to an FFT of our selected block length
 def filtfft(filt, blocklen):
     return sps.freqz(filt[0], filt[1], blocklen, whole=1)[1]