From 280b326947173acddaac3190b5b53feb7fcc48fb Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sun, 21 Jun 2020 11:01:50 -0700 Subject: [PATCH 1/9] a little bit of code movement --- lddecode/core.py | 285 ++++++++++++++++++++-------------------------- lddecode/utils.py | 21 ++++ 2 files changed, 147 insertions(+), 159 deletions(-) diff --git a/lddecode/core.py b/lddecode/core.py index 87fc7d476..6bb86ed1a 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -89,6 +89,9 @@ def calclinelen(SP, mult, mhz): 'hsyncPulseUS': 4.7, 'eqPulseUS': 2.3, 'vsyncPulseUS': 27.1, + + # What 0 IRE/0V should be in digitaloutput + 'outputZero': 1024, } # In color NTSC, the line period was changed from 63.5 to 227.5 color cycles, @@ -134,6 +137,9 @@ def calclinelen(SP, mult, mhz): 'hsyncPulseUS': 4.7, 'eqPulseUS': 2.35, 'vsyncPulseUS': 27.3, + + # What 0 IRE/0V should be in digitaloutput + 'outputZero': 256, } SysParams_PAL['outlinelen'] = calclinelen(SysParams_PAL, 4, 'fsc_mhz') @@ -760,7 +766,7 @@ def audio_phase2(self, field_audio): output_audio2[:tmp.shape[0]] = tmp - end = field_audio.shape[0] #// filterset['audio_fdiv2'] + end = field_audio.shape[0] askip = 512 # length of filters that needs to be chopped out of the ifft sjump = self.blocklen - (askip * self.Filters['audio_fdiv2']) @@ -1117,10 +1123,6 @@ def setparams(self, params): # Apply params to the core thread, so they match up with the decoders self.apply_newparams(params) -@njit -def dsa_rescale(infloat): - return int(np.round(infloat * 32767 / 150000)) - # Downscales to 16bit/44.1khz. It might be nice when analog audio is better to support 24/96, # but if we only support one output type, matching CD audio/digital sound is greatly preferable. def downscale_audio(audio, lineinfo, rf, linecount, timeoffset = 0, freq = 48000.0, scale=64): @@ -1273,6 +1275,15 @@ def usectooutpx(self, x): def outpxtousec(self, x): return x / self.rf.SysParams['outfreq'] + def hz_to_output(self, input): + reduced = (input - self.rf.SysParams['ire0']) / self.rf.SysParams['hz_ire'] + reduced -= self.rf.SysParams['vsync_ire'] + + return np.uint16(np.clip((reduced * self.out_scale) + self.rf.SysParams['outputZero'], 0, 65535) + 0.5) + + def output_to_ire(self, output): + return ((output - self.rf.SysParams['outputZero']) / self.out_scale) + self.rf.SysParams['vsync_ire'] + def lineslice_tbc(self, l, begin = None, length = None, linelocs = None, keepphase = False): ''' return a slice corresponding with pre-TBC line l ''' @@ -2096,15 +2107,6 @@ def calc_burstmedian(self): return nb_median(burstlevel / self.rf.SysParams['hz_ire']) - def hz_to_output(self, input): - reduced = (input - self.rf.SysParams['ire0']) / self.rf.SysParams['hz_ire'] - reduced -= self.rf.SysParams['vsync_ire'] - - return np.uint16(np.clip((reduced * self.out_scale) + 256, 0, 65535) + 0.5) - - def output_to_ire(self, output): - return ((output- 256.0) / self.out_scale) + self.rf.SysParams['vsync_ire'] - def downscale(self, final = False, *args, **kwargs): # For PAL, each field starts with the line containing the first full VSYNC pulse dsout, dsaudio, dsefm = super(FieldPAL, self).downscale(*args, **kwargs) @@ -2145,21 +2147,110 @@ def process(self): #self.downscale(final=True) -# These classes extend Field to do PAL/NTSC specific TBC features. +# ... now for NTSC +# This class is a very basic NTSC 1D color decoder, used for alignment etc +# XXX: make this callable earlier on and/or merge into FieldNTSC? +class CombNTSC: + ''' *partial* NTSC comb filter class - only enough to do VITS calculations ATM ''' + + def __init__(self, field): + self.field = field + self.cbuffer = self.buildCBuffer() -# Hotspots found in profiling are refactored here and boosted by numba's jit. -@njit(cache=True) -def clb_findnextburst(burstarea, i, endburstarea, threshold): - for j in range(i, endburstarea): - if np.abs(burstarea[j]) > threshold: - return j, burstarea[j] + def getlinephase(self, line): + ''' determine if a line has positive color burst phase ''' + fieldID = self.field.fieldPhaseID + + fieldPositivePhase = (fieldID == 1) | (fieldID == 4) + + return fieldPositivePhase if ((line % 2) == 0) else not fieldPositivePhase - return None + def buildCBuffer(self, subset = None): + ''' + prev_field: Compute values for previous field + subset: a slice computed by lineslice_tbc (default: whole field) + + NOTE: first and last two returned values will be zero, so slice accordingly + ''' + + data = self.field.dspicture + + if subset: + data = data[subset] + + # this is a translation of this code from tools/ld-chroma-decoder/comb.cpp: + # + # for (qint32 h = configuration.activeVideoStart; h < configuration.activeVideoEnd; h++) { + # qreal tc1 = (((line[h + 2] + line[h - 2]) / 2) - line[h]); + + fldata = data.astype(np.float32) + cbuffer = np.zeros_like(fldata) + + cbuffer[2:-2] = (fldata[:-4] + fldata[4:]) / 2 + cbuffer[2:-2] -= fldata[2:-2] + + return cbuffer -@njit(cache=True) -def clb_subround(x): - # Yes, this was a hotspot. - return np.round(x) - x + def splitIQ(self, cbuffer, line = 0): + ''' + NOTE: currently? only works on one line + + This returns normalized I and Q arrays, each one half the length of cbuffer + ''' + linephase = self.getlinephase(line) + + sq = cbuffer[::2].copy() + si = cbuffer[1::2].copy() + + if not linephase: + si[0::2] = -si[0::2] + sq[1::2] = -sq[1::2] + else: + si[1::2] = -si[1::2] + sq[0::2] = -sq[0::2] + + return si, sq + + def calcLine19Info(self, comb_field2 = None): + ''' returns color burst phase (ideally 147 degrees) and (unfiltered!) SNR ''' + + # Don't need the whole line here, but start at 0 to definitely have an even # + l19_slice = self.field.lineslice_tbc(19, 0, 40) + l19_slice_i70 = self.field.lineslice_tbc(19, 14, 18) + + # fail out if there is obviously bad data + if not ((np.max(self.field.output_to_ire(self.field.dspicture[l19_slice_i70])) < 100) and + (np.min(self.field.output_to_ire(self.field.dspicture[l19_slice_i70])) > 40)): + #logging.info("WARNING: line 19 data incorrect") + #logging.info(np.max(self.field.output_to_ire(self.field.dspicture[l19_slice_i70])), np.min(self.field.output_to_ire(self.field.dspicture[l19_slice_i70]))) + return None, None, None + + cbuffer = self.cbuffer[l19_slice] + if comb_field2 is not None: + # fail out if there is obviously bad data + if not ((np.max(self.field.output_to_ire(comb_field2.field.dspicture[l19_slice_i70])) < 100) and + (np.min(self.field.output_to_ire(comb_field2.field.dspicture[l19_slice_i70])) > 40)): + return None, None, None + + cbuffer -= comb_field2.cbuffer[l19_slice] + cbuffer /= 2 + + si, sq = self.splitIQ(cbuffer, 19) + + sl = slice(110,230) + cdata = np.sqrt((si[sl] ** 2.0) + (sq[sl] ** 2.0)) + + phase = np.arctan2(np.mean(si[sl]),np.mean(sq[sl]))*180/np.pi + if phase < 0: + phase += 360 + + # compute SNR + signal = np.mean(cdata) + noise = np.std(cdata) + + snr = 20 * np.log10(signal / noise) + + return signal / (2 * self.field.out_scale), phase, snr class FieldNTSC(Field): @@ -2168,7 +2259,6 @@ def get_burstlevel(self, l, linelocs = None): return rms(burstarea) * np.sqrt(2) - def compute_line_bursts(self, linelocs, line): ''' Compute the zero crossing for the given line using calczc @@ -2211,13 +2301,14 @@ def compute_line_bursts(self, linelocs, line): zc_bursts_t = np.zeros(64, dtype=np.float) zc_bursts_n = np.zeros(64, dtype=np.float) + # this subroutine is in utils.py, broken out so it can be JIT'd i = clb_findnextburst(burstarea, 0, len(burstarea) - 1, threshold) zc = 0 while i is not None and zc is not None: zc = calczc(burstarea, i[0], 0) if zc is not None: - zc_burst = clb_subround((bstart+zc-s_rem) / zcburstdiv) + zc_burst = distance_from_round((bstart+zc-s_rem) / zcburstdiv) if i[1] < 0: zc_bursts_t[numpos] = zc_burst numpos = numpos + 1 @@ -2333,15 +2424,6 @@ def refine_linelocs_burst(self, linelocs = None): return linelocs_adj, burstlevel - def hz_to_output(self, input): - reduced = (input - self.rf.SysParams['ire0']) / self.rf.SysParams['hz_ire'] - reduced -= self.rf.SysParams['vsync_ire'] - - return np.uint16(np.clip((reduced * self.out_scale) + 1024, 0, 65535) + 0.5) - - def output_to_ire(self, output): - return ((output - 1024) / self.out_scale) + self.rf.SysParams['vsync_ire'] - def downscale(self, lineoffset = 0, final = False, *args, **kwargs): if final == False: if 'audio' in kwargs: @@ -2397,108 +2479,6 @@ def process(self): self.wowfactor = self.computewow(self.linelocs) -class CombNTSC: - ''' *partial* NTSC comb filter class - only enough to do VITS calculations ATM ''' - - def __init__(self, field): - self.field = field - self.cbuffer = self.buildCBuffer() - - def getlinephase(self, line): - ''' determine if a line has positive color burst phase ''' - fieldID = self.field.fieldPhaseID - - fieldPositivePhase = (fieldID == 1) | (fieldID == 4) - - return fieldPositivePhase if ((line % 2) == 0) else not fieldPositivePhase - - def buildCBuffer(self, subset = None): - ''' - prev_field: Compute values for previous field - subset: a slice computed by lineslice_tbc (default: whole field) - - NOTE: first and last two returned values will be zero, so slice accordingly - ''' - - data = self.field.dspicture - - if subset: - data = data[subset] - - # this is a translation of this code from tools/ld-chroma-decoder/comb.cpp: - # - # for (qint32 h = configuration.activeVideoStart; h < configuration.activeVideoEnd; h++) { - # qreal tc1 = (((line[h + 2] + line[h - 2]) / 2) - line[h]); - - fldata = data.astype(np.float32) - cbuffer = np.zeros_like(fldata) - - cbuffer[2:-2] = (fldata[:-4] + fldata[4:]) / 2 - cbuffer[2:-2] -= fldata[2:-2] - - return cbuffer - - def splitIQ(self, cbuffer, line = 0): - ''' - NOTE: currently? only works on one line - - This returns normalized I and Q arrays, each one half the length of cbuffer - ''' - linephase = self.getlinephase(line) - - sq = cbuffer[::2].copy() - si = cbuffer[1::2].copy() - - if not linephase: - si[0::2] = -si[0::2] - sq[1::2] = -sq[1::2] - else: - si[1::2] = -si[1::2] - sq[0::2] = -sq[0::2] - - return si, sq - - def calcLine19Info(self, comb_field2 = None): - ''' returns color burst phase (ideally 147 degrees) and (unfiltered!) SNR ''' - - # Don't need the whole line here, but start at 0 to definitely have an even # - l19_slice = self.field.lineslice_tbc(19, 0, 40) - l19_slice_i70 = self.field.lineslice_tbc(19, 14, 18) - - # fail out if there is obviously bad data - if not ((np.max(self.field.output_to_ire(self.field.dspicture[l19_slice_i70])) < 100) and - (np.min(self.field.output_to_ire(self.field.dspicture[l19_slice_i70])) > 40)): - #logging.info("WARNING: line 19 data incorrect") - #logging.info(np.max(self.field.output_to_ire(self.field.dspicture[l19_slice_i70])), np.min(self.field.output_to_ire(self.field.dspicture[l19_slice_i70]))) - return None, None, None - - cbuffer = self.cbuffer[l19_slice] - if comb_field2 is not None: - # fail out if there is obviously bad data - if not ((np.max(self.field.output_to_ire(comb_field2.field.dspicture[l19_slice_i70])) < 100) and - (np.min(self.field.output_to_ire(comb_field2.field.dspicture[l19_slice_i70])) > 40)): - return None, None, None - - cbuffer -= comb_field2.cbuffer[l19_slice] - cbuffer /= 2 - - si, sq = self.splitIQ(cbuffer, 19) - - sl = slice(110,230) - cdata = np.sqrt((si[sl] ** 2.0) + (sq[sl] ** 2.0)) - - phase = np.arctan2(np.mean(si[sl]),np.mean(sq[sl]))*180/np.pi - if phase < 0: - phase += 360 - - # compute SNR - signal = np.mean(cdata) - noise = np.std(cdata) - - snr = 20 * np.log10(signal / noise) - - return signal / (2 * self.field.out_scale), phase, snr - class LDdecode: def __init__(self, fname_in, fname_out, freader, analog_audio = 0, digital_audio = False, system = 'NTSC', doDOD = True, threads=4, extra_options = {}): @@ -2614,29 +2594,17 @@ def roughseek(self, location, isField = True): else: self.fdoffset = location - def checkMTF_calc(self, field): - if not self.isCLV and self.frameNumber is not None: - newmtf = 1 - (self.frameNumber / 10000) if self.frameNumber < 10000 else 0 - oldmtf = self.mtf_level - self.mtf_level = newmtf - - if np.abs(newmtf - oldmtf) > .1: # redo field if too much adjustment - #logging.info(newmtf, oldmtf, field.cavFrame) - return False - - return True - def checkMTF(self, field, pfield = None): - if not self.autoMTF: - return self.checkMTF_calc(field) - oldmtf = self.mtf_level - if (len(self.bw_ratios) == 0): - return True + if not self.autoMTF: + self.mtf_level = 1 - (self.frameNumber / 10000) if self.frameNumber < 10000 else 0 + else: + if (len(self.bw_ratios) == 0): + return True - # scale for NTSC - 1.1 to 1.55 - self.mtf_level = np.clip((np.mean(self.bw_ratios) - 1.08) / .38, 0, 1) + # scale for NTSC - 1.1 to 1.55 + self.mtf_level = np.clip((np.mean(self.bw_ratios) - 1.08) / .38, 0, 1) return np.abs(self.mtf_level - oldmtf) < .05 @@ -2994,7 +2962,6 @@ def computeMetrics(self, f, fp = None, verbose = False): metrics['bPSNR'] = self.calcpsnr(f, bl_slicetbc) if 'whiteRFLevel' in metrics: - #metrics['syncToWhiteRFRatio'] = metrics['syncRFLevel'] / metrics['whiteRFLevel'] metrics['blackToWhiteRFRatio'] = metrics['blackLineRFLevel'] / metrics['whiteRFLevel'] outputkeys = metrics.keys() if verbose else ['wSNR', 'bPSNR'] diff --git a/lddecode/utils.py b/lddecode/utils.py index ae87db8e1..0547aa313 100644 --- a/lddecode/utils.py +++ b/lddecode/utils.py @@ -805,6 +805,27 @@ def db_to_lev(db): def lev_to_db(rlev): return 20 * np.log10(rlev) +# moved from core.py +@njit +def dsa_rescale(infloat): + return int(np.round(infloat * 32767 / 150000)) + +# Hotspot subroutines in FieldNTSC's compute_line_bursts function, +# removed so that they can be JIT'd + +@njit(cache=True) +def clb_findnextburst(burstarea, i, endburstarea, threshold): + for j in range(i, endburstarea): + if np.abs(burstarea[j]) > threshold: + return j, burstarea[j] + + return None + +@njit(cache=True) +def distance_from_round(x): + # Yes, this was a hotspot. + return np.round(x) - x + # Write the .tbc.json file (used by lddecode and notebooks) def write_json(ldd, outname): jsondict = ldd.build_json(ldd.curfield) From 314d5b59892ad4afda5bc6327af27cd26e0ea9b0 Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sun, 21 Jun 2020 11:38:45 -0700 Subject: [PATCH 2/9] a little more pruning --- lddecode/core.py | 54 ++++++++++++++++++--------------------- lddecode/plot_utils.py | 56 ++++++++++++++++++++++++++++++++++++++++ lddecode/utils.py | 58 ------------------------------------------ 3 files changed, 81 insertions(+), 87 deletions(-) diff --git a/lddecode/core.py b/lddecode/core.py index 6bb86ed1a..48b01c5cd 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -323,7 +323,7 @@ def __init__(self, inputfreq = 40, system = 'NTSC', blocklen = 32*1024, decode_d self.freq = freq self.freq_half = freq / 2 self.freq_hz = self.freq * 1000000 - self.freq_hz_half = self.freq * 1000000 / 2 + self.freq_hz_half = self.freq_hz / 2 self.mtf_mult = 1.0 self.mtf_offset = 0 @@ -360,7 +360,8 @@ def __init__(self, inputfreq = 40, system = 'NTSC', blocklen = 32*1024, decode_d self.computefilters() - # The 0.5mhz filter is rolled back, so there are a few unusable samples at the end + # The 0.5mhz filter is rolled back to align with the data, so there + # are a few unusable samples at the end. self.blockcut_end = self.Filters['F05_offset'] def computefilters(self): @@ -504,6 +505,10 @@ def audio_fdslice(self, freqdomain): def audio_fdslice2(self, freqdomain): return np.concatenate([freqdomain[self.Filters['audio_fdslice2_lo']], freqdomain[self.Filters['audio_fdslice2_hi']]]) + def compute_deemp_audio2(self, dfreq): + adeemp_b, adeemp_a = sps.butter(1, [(1000000/dfreq)/(self.Filters['freq_aud2']/2)], btype='lowpass') + return filtfft([adeemp_b, adeemp_a], self.blocklen // self.Filters['audio_fdiv2']) + def computeaudiofilters(self): SF = self.Filters SP = self.SysParams @@ -556,17 +561,10 @@ def computeaudiofilters(self): # XXX: This probably needs further tuning, but more or less flattens the 20hz-20khz response # on both PAL and NTSC - d75freq = 1000000/(2*pi*62) - adeemp_b, adeemp_a = sps.butter(1, [d75freq/(SF['freq_aud2']/2)], btype='lowpass') - addemp2lp = filtfft([adeemp_b, adeemp_a], self.blocklen // audio_fdiv2) - - dxfreq = 1000000/(2*pi*45) - adeemp_b, adeemp_a = sps.butter(1, [dxfreq/(SF['freq_aud2']/2)], btype='highpass') - addemp2hp1 = filtfft([adeemp_b, adeemp_a], self.blocklen // audio_fdiv2) - - dxfreq = 1000000/(2*pi*8) - adeemp_b, adeemp_a = sps.butter(1, [dxfreq/(SF['freq_aud2']/2)], btype='highpass') - addemp2hp2 = filtfft([adeemp_b, adeemp_a], self.blocklen // audio_fdiv2) + # (and no, I don't know where those frequencies come from) + addemp2lp = self.compute_deemp_audio2(2*pi*62) # 2567hz + addemp2hp1 = self.compute_deemp_audio2(2*pi*45) # 3536hz + addemp2hp2 = self.compute_deemp_audio2(2*pi*8) # 19894hz (i.e. cutoff?) SF['audio_deemp2'] = addemp2lp + (addemp2hp1 * .14) + (addemp2hp2 * .29) @@ -2507,24 +2505,23 @@ def __init__(self, fname_in, fname_out, freader, analog_audio = 0, digital_audio if fname_out is not None: self.outfile_video = open(fname_out + '.tbc', 'wb') - #self.outfile_json = open(fname_out + '.json', 'wb') self.outfile_audio = open(fname_out + '.pcm', 'wb') if self.analog_audio else None self.outfile_efm = open(fname_out + '.efm', 'wb') if self.digital_audio else None + + if digital_audio: + # feed EFM stream into ld-ldstoefm + self.efm_pll = efm_pll.EFM_PLL() else: self.outfile_video = None self.outfile_audio = None self.outfile_efm = None - if fname_out is not None and digital_audio: - # feed EFM stream into ld-ldstoefm - self.efm_pll = efm_pll.EFM_PLL() - self.fname_out = fname_out self.firstfield = None # In frame output mode, the first field goes here self.fieldloc = 0 - self.internalerrors = 0 # exceptions seen + self.internalerrors = [] # exceptions seen # if option is missing, get returns None @@ -2580,7 +2577,7 @@ def __del__(self): def close(self): ''' deletes all open files, so it's possible to pickle an LDDecode object ''' - # use setattr to force file closure + # use setattr to force file closure by unlinking the objects for outfiles in ['infile', 'outfile_video', 'outfile_audio', 'outfile_json', 'outfile_efm']: setattr(self, outfiles, None) @@ -2589,10 +2586,9 @@ def close(self): def roughseek(self, location, isField = True): self.prevPhaseID = None + self.fdoffset = location if isField: - self.fdoffset = location * self.bytes_per_field - else: - self.fdoffset = location + self.fdoffset *= self.bytes_per_field def checkMTF(self, field, pfield = None): oldmtf = self.mtf_level @@ -2653,8 +2649,6 @@ def writeout(self, dataset): if self.digital_audio == True: efm_out = self.efm_pll.process(efm) self.outfile_efm.write(efm_out.tobytes()) - else: - efm = None def decodefield(self, initphase = False): ''' returns field object if valid, and the offset to the next decode ''' @@ -2816,7 +2810,6 @@ def decodeBCD(bcd): try: self.clvMinutes = decodeBCD(l & 0xff) + (decodeBCD((l >> 16) & 0xf) * 60) self.isCLV = True - #logging.info('CLV', mins) except ValueError: pass elif (l & 0xf00000) == 0xf00000: # CAV frame @@ -2839,13 +2832,16 @@ def decodeBCD(bcd): pass if self.clvMinutes is not None: + minute_seconds = self.clvMinutes * 60 + if self.clvSeconds is not None: # newer CLV - return (((self.clvMinutes * 60) + self.clvSeconds) * self.clvfps) + self.clvFrameNum + # XXX: does not auto-decrement for skip frames + return ((minute_seconds + self.clvSeconds) * self.clvfps) + self.clvFrameNum else: self.earlyCLV = True - return (self.clvMinutes * 60) + return minute_seconds - return None #seeking won't work w/minutes only + return None # seeking won't work w/minutes only def calcsnr(self, f, snrslice): data = f.output_to_ire(f.dspicture[snrslice]) diff --git a/lddecode/plot_utils.py b/lddecode/plot_utils.py index edc0111bc..68fea75c0 100644 --- a/lddecode/plot_utils.py +++ b/lddecode/plot_utils.py @@ -10,6 +10,7 @@ import scipy as sp import scipy.signal as sps +import matplotlib import matplotlib.pyplot as plt # To support image displays @@ -17,6 +18,61 @@ import IPython.display from IPython.display import HTML +def todb(y, zero = False): + db = 20 * np.log10(np.abs(y)) + if zero: + return db - np.max(db) + else: + return db + +def plotfilter_wh(w, h, freq, zero_base = False): + db = todb(h, zero_base) + + above_m3 = None + for i in range(1, len(w)): + if (db[i] >= -10) and (db[i - 1] < -10): + print(">-10db crossing at ", w[i]) + if (db[i] >= -3) and (db[i - 1] < -3): + print(">-3db crossing at ", w[i]) + above_m3 = i + if (db[i] < -3) and (db[i - 1] >= -3): + if above_m3 is not None: + peak_index = np.argmax(db[above_m3:i]) + above_m3 + print("peak at ", w[peak_index], db[peak_index]) + print("<-3db crossing at ", w[i]) + if (db[i] >= 3) and (db[i - 1] < 3): + print(">3db crossing at ", w[i]) + + fig, ax1 = plt.subplots(1, 1, sharex=True) + ax1.set_title('Digital filter frequency response') + + ax1.plot(w, db, 'b') + ax1.set_ylabel('Amplitude [dB]', color='b') + ax1.set_xlabel('Frequency [rad/sample]') + + ax2 = ax1.twinx() + angles = np.unwrap(np.angle(h)) + ax2.plot(w, angles, 'g') + ax2.set_ylabel('Angle (radians)', color='g') + + plt.grid() + plt.axis('tight') + plt.show() + + return None + +def plotfilter(B, A, dfreq = None, freq = 40, zero_base = False): + if dfreq is None: + dfreq = freq / 2 + + w, h = sps.freqz(B, A, whole=True, worN=4096) + w = np.arange(0, freq, freq / len(h)) + + keep = int((dfreq / freq) * len(h)) + + return plotfilter_wh(w[1:keep], h[1:keep], freq, zero_base) + + pi = np.pi tau = np.pi * 2 diff --git a/lddecode/utils.py b/lddecode/utils.py index 0547aa313..ff611f07b 100644 --- a/lddecode/utils.py +++ b/lddecode/utils.py @@ -19,64 +19,6 @@ import scipy as sp import scipy.signal as sps -# plotting -import matplotlib -import matplotlib.pyplot as plt - -def todb(y, zero = False): - db = 20 * np.log10(np.abs(y)) - if zero: - return db - np.max(db) - else: - return db - -def plotfilter_wh(w, h, freq, zero_base = False): - db = todb(h, zero_base) - - above_m3 = None - for i in range(1, len(w)): - if (db[i] >= -10) and (db[i - 1] < -10): - print(">-10db crossing at ", w[i]) - if (db[i] >= -3) and (db[i - 1] < -3): - print(">-3db crossing at ", w[i]) - above_m3 = i - if (db[i] < -3) and (db[i - 1] >= -3): - if above_m3 is not None: - peak_index = np.argmax(db[above_m3:i]) + above_m3 - print("peak at ", w[peak_index], db[peak_index]) - print("<-3db crossing at ", w[i]) - if (db[i] >= 3) and (db[i - 1] < 3): - print(">3db crossing at ", w[i]) - - fig, ax1 = plt.subplots(1, 1, sharex=True) - ax1.set_title('Digital filter frequency response') - - ax1.plot(w, db, 'b') - ax1.set_ylabel('Amplitude [dB]', color='b') - ax1.set_xlabel('Frequency [rad/sample]') - - ax2 = ax1.twinx() - angles = np.unwrap(np.angle(h)) - ax2.plot(w, angles, 'g') - ax2.set_ylabel('Angle (radians)', color='g') - - plt.grid() - plt.axis('tight') - plt.show() - - return None - -def plotfilter(B, A, dfreq = None, freq = 40, zero_base = False): - if dfreq is None: - dfreq = freq / 2 - - w, h = sps.freqz(B, A, whole=True, worN=4096) - w = np.arange(0, freq, freq / len(h)) - - keep = int((dfreq / freq) * len(h)) - - return plotfilter_wh(w[1:keep], h[1:keep], freq, zero_base) - from scipy import interpolate # This uses numpy's interpolator, which works well enough From bd3be3f6d5181eb6808d3ff8fcb9517df9d4237b Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sun, 21 Jun 2020 12:06:47 -0700 Subject: [PATCH 3/9] yet more pruning --- lddecode/core.py | 76 ++++++++++++++++-------------------------------- 1 file changed, 25 insertions(+), 51 deletions(-) diff --git a/lddecode/core.py b/lddecode/core.py index 48b01c5cd..cb9fa302c 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -370,16 +370,12 @@ def computefilters(self): self.computevideofilters() # This is > 0 because decode_analog_audio is in khz. - if self.decode_analog_audio > 0: + if self.decode_analog_audio != 0: self.computeaudiofilters() if self.decode_digital_audio: self.computeefmfilter() - # This high pass filter is intended to detect RF dropouts - Frfhpf = sps.butter(1, [10/self.freq_half], btype='highpass') - self.Filters['Frfhpf'] = filtfft(Frfhpf, self.blocklen) - self.computedelays() def computeefmfilter(self): @@ -409,7 +405,6 @@ def computeefmfilter(self): p_interp = spi.interp1d(freqs, phase, kind="cubic") nonzero_bins = int(freqs[-1] / freq_per_bin) + 1 - #print(nonzero_bins) bin_freqs = np.arange(nonzero_bins) * freq_per_bin bin_amp = a_interp(bin_freqs) @@ -418,9 +413,6 @@ def computeefmfilter(self): # Scale by the amplitude, rotate by the phase coeffs[:nonzero_bins] = bin_amp * (np.cos(bin_phase) + (complex(0, -1) * np.sin(bin_phase))) - # TODO: figure out this flip, the filter is asymmetric but works better that way. huh. - #coeffs[-nonzero_bins:] = np.flip(coeffs[:nonzero_bins], 0) - self.Filters['Fefm'] = coeffs * 8 def computevideofilters(self): @@ -431,6 +423,10 @@ def computevideofilters(self): SP = self.SysParams DP = self.DecoderParams + # This high pass filter is intended to detect RF dropouts + Frfhpf = sps.butter(1, [10/self.freq_half], btype='highpass') + self.Filters['Frfhpf'] = filtfft(Frfhpf, self.blocklen) + # First phase FFT filtering # MTF filter section @@ -482,7 +478,6 @@ def computevideofilters(self): # Post processing: lowpass filter + deemp SF['FVideo'] = SF['Fvideo_lpf'] * SF['Fdeemp'] - #SF['FVideo'] = SF['Fdeemp'] # additional filters: 0.5mhz and color burst # Using an FIR filter here to get a known delay @@ -843,23 +838,18 @@ def computedelays(self, mtf_level = 0): # XXX: sync detector does NOT reflect actual sync detection, just regular filtering @ sync level # (but only regular filtering is needed for DOD) - dgap_sync = calczc(fakedecode['video']['demod'], 1500, rf.iretohz(rf.SysParams['vsync_ire'] / 2), count=512) - 1500 - dgap_white = calczc(fakedecode['video']['demod'], 3000, rf.iretohz(50), count=512) - 3000 - dgap_rot = calczc(fakedecode['video']['demod'], 6000, rf.iretohz(-10), count=512) - 6000 - rf.delays = {} - # factor in the 1k or so block cut as well, since we never *just* do demodblock - rf.delays['video_sync'] = dgap_sync #- self.blockcut - rf.delays['video_white'] = dgap_white #- self.blockcut - rf.delays['video_rot'] = int(np.round(dgap_rot)) #- self.blockcut - + rf.delays['video_sync'] = calczc(fakedecode['video']['demod'], 1500, rf.iretohz(rf.SysParams['vsync_ire'] / 2), count=512) - 1500 + rf.delays['video_white'] = calczc(fakedecode['video']['demod'], 3000, rf.iretohz(50), count=512) - 3000 + rf.delays['video_rot'] = int(np.round(calczc(fakedecode['video']['demod'], 6000, rf.iretohz(-10), count=512) - 6000)) + fdec_raw = fakedecode['video']['demod_raw'] rf.limits = {} rf.limits['sync'] = (np.min(fdec_raw[1400:2800]), np.max(fdec_raw[1400:2800])) rf.limits['viewable'] = (np.min(fdec_raw[2900:6000]), np.max(fdec_raw[2900:6000])) - return fakedecode, dgap_sync, dgap_white + return fakedecode class DemodCache: def __init__(self, rf, infile, loader, cachesize = 256, num_worker_threads=6, MTF_tolerance = .05): @@ -2334,13 +2324,11 @@ def compute_burst_offsets(self, linelocs): if (len(zc_bursts[l][True]) == 0) or (len(zc_bursts[l][False]) == 0): continue + + even_line = not (l % 2) - if not (l % 2): - bursts['even'].append(zc_bursts[l][True]) - bursts['odd'].append(zc_bursts[l][False]) - else: - bursts['even'].append(zc_bursts[l][False]) - bursts['odd'].append(zc_bursts[l][True]) + bursts['even'].append(zc_bursts[l][True if even_line else False]) + bursts['odd'].append(zc_bursts[l][False if even_line else True]) bursts_arr = {} bursts_arr[True] = np.concatenate(bursts['even']) @@ -2718,10 +2706,7 @@ def readfield(self, initphase = False): #logging.info(metrics['blackToWhiteRFRatio'], np.mean(self.bw_ratios)) - redo = False - - if not self.checkMTF(f, self.prevfield): - redo = True + redo = not self.checkMTF(f, self.prevfield) # Perform AGC changes on first fields only to prevent luma mismatch intra-field if self.useAGC and f.isFirstField: @@ -2786,24 +2771,17 @@ def decodeBCD(bcd): raise ValueError("Non-decimal BCD digit") return (10 * decodeBCD(bcd >> 4)) + digit - # With newer versions of the duplicator software, there may - # only be one *field* with leadout, so handle each field - # in turn. (issue #414) - for i in [f1.linecode, f2.linecode]: - leadoutCount = 0 - for v in i: - if v is not None and v == 0x80eeee: - leadoutCount += 1 - - if leadoutCount == 2: - self.leadOut = True + leadoutCount = 0 for l in f1.linecode + f2.linecode: if l is None: continue if l == 0x80eeee: # lead-out reached - pass + leadoutCount += 1 + # Require two leadouts, since there may only be one field in the raw data w/it + if leadoutCount == 2: + self.leadOut = True elif l == 0x88ffff: # lead-in self.leadIn = True elif (l & 0xf0dd00) == 0xf0dd00: # CLV minutes/hours @@ -2843,21 +2821,17 @@ def decodeBCD(bcd): return None # seeking won't work w/minutes only - def calcsnr(self, f, snrslice): - data = f.output_to_ire(f.dspicture[snrslice]) + def calcsnr(self, f, snrslice, psnr = False): + # if dspicture isn't converted to float, this underflows at -40IRE + data = f.output_to_ire(f.dspicture[snrslice].astype(float)) - signal = np.mean(data) + signal = np.mean(data) if not psnr else 100 noise = np.std(data) return 20 * np.log10(signal / noise) def calcpsnr(self, f, snrslice): - # if dspicture isn't converted to float, this underflows at -40IRE - data = f.output_to_ire(f.dspicture[snrslice].astype(float)) - - noise = np.std(data) - - return 20 * np.log10(100 / noise) + return self.calcsnr(f, snrslice, psnr=True) def computeMetricsPAL(self, metrics, f, fp = None): From 163f750b9be9e78b6d04d8833c6bd4f681180fa3 Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sat, 27 Jun 2020 21:44:46 -0700 Subject: [PATCH 4/9] fix vroopdark --- lddecode/core.py | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/lddecode/core.py b/lddecode/core.py index cb9fa302c..18e5258cd 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -1767,10 +1767,12 @@ def refine_linelocs_hsync(self): linelocs2[i] = zc2 else: self.linebad[i] = True - linelocs2[i] = self.linelocs1[i] # don't use the computed value here if it's bad else: self.linebad[i] = True + if self.linebad[i]: + linelocs2[i] = self.linelocs1[i] # don't use the computed value here if it's bad + return linelocs2 def compute_deriv_error(self, linelocs, baserr): @@ -1795,6 +1797,13 @@ def fix_badlines(self, linelocs_in, linelocs_backup_in = None): badlines = np.isnan(linelocs) linelocs[badlines] = linelocs_backup[badlines] + # compute mean adjustment for linelocs1 + adjs = [] + for l in range(20, self.linecount - 1): + if l not in self.linebad: + adjs.append(linelocs_in[l] - self.linelocs1[l]) + adjmean = np.mean(adjs) + for l in np.where(self.linebad)[0]: prevgood = l - 1 nextgood = l + 1 @@ -1805,12 +1814,14 @@ def fix_badlines(self, linelocs_in, linelocs_backup_in = None): while nextgood < len(linelocs) and self.linebad[nextgood]: nextgood += 1 - #print(l, prevgood, nextgood) - - if prevgood >= 0 and nextgood < len(linelocs): + if prevgood > 0 and nextgood <= self.linecount: gap = (linelocs[nextgood] - linelocs[prevgood]) / (nextgood - prevgood) linelocs[l] = (gap * (l - prevgood)) + linelocs[prevgood] - + else: + linelocs[l] = self.linelocs1[l] + adjmean + + #print(l, prevgood, nextgood, self.linelocs1[l] - self.linelocs1[l - 1], linelocs[l] - linelocs[l-1]) + return linelocs def computewow(self, lineinfo): @@ -2617,7 +2628,7 @@ def detectLevels(self, field): ire = field.rf.hztoire(hz) if ire < field.rf.SysParams['vsync_ire'] / 2: sync_hzs.append(hz) - else: + elif ire < 10: ire0_hzs.append(hz) return np.mean(sync_hzs), np.mean(ire0_hzs) @@ -2713,6 +2724,8 @@ def readfield(self, initphase = False): sync_hz, ire0_hz = self.detectLevels(f) sync_ire_diff = np.abs(self.rf.hztoire(sync_hz) - self.rf.SysParams['vsync_ire']) + #print(sync_hz, ire0_hz, sync_ire_diff) + if (sync_ire_diff > 2) or (np.abs(self.rf.hztoire(ire0_hz)) > 2): redo = True From 9191bc8a0f0a1a1c04f08e7e8e4d9f70cbf3e220 Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sat, 27 Jun 2020 21:57:54 -0700 Subject: [PATCH 5/9] fix glitchbright (and would fix vroopdark too) --- lddecode/core.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/lddecode/core.py b/lddecode/core.py index 18e5258cd..79dff6033 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -2609,7 +2609,7 @@ def detectLevels(self, field): # Build a list of each half-line's average hlevels = [] - for l in range(1,8): + for l in range(2,8): lsa = field.lineslice(l, 10, 10) lsb = field.lineslice(l, 40, 10) @@ -2619,6 +2619,8 @@ def detectLevels(self, field): hlevels.append(np.median(field.data['video']['demod_05'][lsa]) / adj) hlevels.append(np.median(field.data['video']['demod_05'][lsb]) / adj) + + #print(hlevels[-2:]) # Now group them by level (either sync or ire 0) and return the means of those sync_hzs = [] @@ -2726,7 +2728,9 @@ def readfield(self, initphase = False): #print(sync_hz, ire0_hz, sync_ire_diff) - if (sync_ire_diff > 2) or (np.abs(self.rf.hztoire(ire0_hz)) > 2): + diff_to_redo = 2 if self.fields_written else .5 + + if (sync_ire_diff > diff_to_redo) or (np.abs(self.rf.hztoire(ire0_hz)) > diff_to_redo): redo = True self.rf.SysParams['ire0'] = ire0_hz From f43ac62fe079300744f060065cca7e14348e7ecf Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sat, 27 Jun 2020 22:00:21 -0700 Subject: [PATCH 6/9] back out vroopdark change, broke pal --- lddecode/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lddecode/core.py b/lddecode/core.py index 79dff6033..733f8b846 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -2609,7 +2609,7 @@ def detectLevels(self, field): # Build a list of each half-line's average hlevels = [] - for l in range(2,8): + for l in range(1,8): lsa = field.lineslice(l, 10, 10) lsb = field.lineslice(l, 40, 10) From 0245ebd51491ce84ccc8bd177dfe5f508a7469bb Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sat, 27 Jun 2020 22:00:39 -0700 Subject: [PATCH 7/9] back out vroopdark change, broke pal --- lddecode/core.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lddecode/core.py b/lddecode/core.py index 733f8b846..164ebba46 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -2609,7 +2609,7 @@ def detectLevels(self, field): # Build a list of each half-line's average hlevels = [] - for l in range(1,8): + for l in range(2,8): lsa = field.lineslice(l, 10, 10) lsb = field.lineslice(l, 40, 10) @@ -2630,7 +2630,7 @@ def detectLevels(self, field): ire = field.rf.hztoire(hz) if ire < field.rf.SysParams['vsync_ire'] / 2: sync_hzs.append(hz) - elif ire < 10: + else: ire0_hzs.append(hz) return np.mean(sync_hzs), np.mean(ire0_hzs) From 959156dc0651d493876fac24da47bd8d643b6fd4 Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sun, 28 Jun 2020 07:04:43 -0700 Subject: [PATCH 8/9] back out last couple of changes --- lddecode/core.py | 31 +++++++------------------------ 1 file changed, 7 insertions(+), 24 deletions(-) diff --git a/lddecode/core.py b/lddecode/core.py index 164ebba46..cb9fa302c 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -1767,12 +1767,10 @@ def refine_linelocs_hsync(self): linelocs2[i] = zc2 else: self.linebad[i] = True + linelocs2[i] = self.linelocs1[i] # don't use the computed value here if it's bad else: self.linebad[i] = True - if self.linebad[i]: - linelocs2[i] = self.linelocs1[i] # don't use the computed value here if it's bad - return linelocs2 def compute_deriv_error(self, linelocs, baserr): @@ -1797,13 +1795,6 @@ def fix_badlines(self, linelocs_in, linelocs_backup_in = None): badlines = np.isnan(linelocs) linelocs[badlines] = linelocs_backup[badlines] - # compute mean adjustment for linelocs1 - adjs = [] - for l in range(20, self.linecount - 1): - if l not in self.linebad: - adjs.append(linelocs_in[l] - self.linelocs1[l]) - adjmean = np.mean(adjs) - for l in np.where(self.linebad)[0]: prevgood = l - 1 nextgood = l + 1 @@ -1814,14 +1805,12 @@ def fix_badlines(self, linelocs_in, linelocs_backup_in = None): while nextgood < len(linelocs) and self.linebad[nextgood]: nextgood += 1 - if prevgood > 0 and nextgood <= self.linecount: + #print(l, prevgood, nextgood) + + if prevgood >= 0 and nextgood < len(linelocs): gap = (linelocs[nextgood] - linelocs[prevgood]) / (nextgood - prevgood) linelocs[l] = (gap * (l - prevgood)) + linelocs[prevgood] - else: - linelocs[l] = self.linelocs1[l] + adjmean - - #print(l, prevgood, nextgood, self.linelocs1[l] - self.linelocs1[l - 1], linelocs[l] - linelocs[l-1]) - + return linelocs def computewow(self, lineinfo): @@ -2609,7 +2598,7 @@ def detectLevels(self, field): # Build a list of each half-line's average hlevels = [] - for l in range(2,8): + for l in range(1,8): lsa = field.lineslice(l, 10, 10) lsb = field.lineslice(l, 40, 10) @@ -2619,8 +2608,6 @@ def detectLevels(self, field): hlevels.append(np.median(field.data['video']['demod_05'][lsa]) / adj) hlevels.append(np.median(field.data['video']['demod_05'][lsb]) / adj) - - #print(hlevels[-2:]) # Now group them by level (either sync or ire 0) and return the means of those sync_hzs = [] @@ -2726,11 +2713,7 @@ def readfield(self, initphase = False): sync_hz, ire0_hz = self.detectLevels(f) sync_ire_diff = np.abs(self.rf.hztoire(sync_hz) - self.rf.SysParams['vsync_ire']) - #print(sync_hz, ire0_hz, sync_ire_diff) - - diff_to_redo = 2 if self.fields_written else .5 - - if (sync_ire_diff > diff_to_redo) or (np.abs(self.rf.hztoire(ire0_hz)) > diff_to_redo): + if (sync_ire_diff > 2) or (np.abs(self.rf.hztoire(ire0_hz)) > 2): redo = True self.rf.SysParams['ire0'] = ire0_hz From a7c970c9d09478f8aa25568ec1de171a2eca5883 Mon Sep 17 00:00:00 2001 From: Chad Page Date: Sun, 28 Jun 2020 07:19:41 -0700 Subject: [PATCH 9/9] redo fixes, pre-test pal this time --- lddecode/core.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/lddecode/core.py b/lddecode/core.py index cb9fa302c..482da189b 100644 --- a/lddecode/core.py +++ b/lddecode/core.py @@ -1767,10 +1767,12 @@ def refine_linelocs_hsync(self): linelocs2[i] = zc2 else: self.linebad[i] = True - linelocs2[i] = self.linelocs1[i] # don't use the computed value here if it's bad else: self.linebad[i] = True + if self.linebad[i]: + linelocs2[i] = self.linelocs1[i] # don't use the computed value here if it's bad + return linelocs2 def compute_deriv_error(self, linelocs, baserr): @@ -1807,7 +1809,8 @@ def fix_badlines(self, linelocs_in, linelocs_backup_in = None): #print(l, prevgood, nextgood) - if prevgood >= 0 and nextgood < len(linelocs): + firstcheck = 0 if self.rf.system == 'PAL' else 1 + if prevgood >= firstcheck and nextgood < (len(linelocs) + self.lineoffset): gap = (linelocs[nextgood] - linelocs[prevgood]) / (nextgood - prevgood) linelocs[l] = (gap * (l - prevgood)) + linelocs[prevgood] @@ -2598,17 +2601,19 @@ def detectLevels(self, field): # Build a list of each half-line's average hlevels = [] - for l in range(1,8): + for l in range(2,8): lsa = field.lineslice(l, 10, 10) lsb = field.lineslice(l, 40, 10) # compute wow adjustment thislinelen = field.linelocs[l + field.lineoffset] - field.linelocs[l + field.lineoffset - 1] adj = field.rf.linelen / thislinelen - + hlevels.append(np.median(field.data['video']['demod_05'][lsa]) / adj) hlevels.append(np.median(field.data['video']['demod_05'][lsb]) / adj) + #print(l, thislinelen, adj, hlevels[-2:]) + # Now group them by level (either sync or ire 0) and return the means of those sync_hzs = [] ire0_hzs = [] @@ -2713,7 +2718,11 @@ def readfield(self, initphase = False): sync_hz, ire0_hz = self.detectLevels(f) sync_ire_diff = np.abs(self.rf.hztoire(sync_hz) - self.rf.SysParams['vsync_ire']) - if (sync_ire_diff > 2) or (np.abs(self.rf.hztoire(ire0_hz)) > 2): + #print(ire0_hz, sync_ire_diff) + + acceptable_diff = 2 if self.fields_written else .5 + + if (sync_ire_diff > acceptable_diff) or (np.abs(self.rf.hztoire(ire0_hz)) > acceptable_diff): redo = True self.rf.SysParams['ire0'] = ire0_hz