UCBerkeleySETI / hyperseti

A SETI / technosignature search code to find intelligent life beyond Earth
https://hyperseti.readthedocs.io
10 stars 4 forks source link

New idea for dedoppler search algorithm: iterative blanking #38

Closed telegraphic closed 2 years ago

telegraphic commented 2 years ago

Works quite well in simple tests! Idea is to find maximum in freq-doppler space, blank signal in freq-time space, then rerun doppler and find next maximum in freq-doppler space. Code gist:

def find_max_idx(data):
    idx = np.unravel_index(np.argmax(data), data.shape)
    v = data[idx]
    return idx, v

def blank_hit(data, metadata, f0, drate, padding=1):
    n_time, n_pol, n_chans = data.shape
    i0     = int((f0 - metadata['frequency_start']) / metadata['frequency_step'])
    i_step =  metadata['time_step'] * drate / metadata['frequency_step']
    i_off  = (i_step * np.arange(n_time)).astype('int64') + i0
    i_time_axis = list(range(n_time))
    for p_off in range(padding):
        data[i_time_axis, :, i_off] = 0
        data[i_time_axis, :, i_off - p_off] = 0
        data[i_time_axis, :, i_off + p_off] = 0
    return data

def find_hits_recursive(data, max_hits=10, thr=19):
    data = np.copy(np.expand_dims(data, axis=1))
    data = normalize(data.astype('float32'), return_space='cpu')

    hits = []
    hit_id = 0
    while hit_id < max_hits:
        dd, _md = dedoppler(data, metadata, max_dd=1.0, return_space='cpu')
        idx, val = find_max_idx(dd.data)
        if val > thr:
            hits.append((*idx, val))
            f_idx  = dd.dims.index('frequency')
            f0     = dd.frequency[idx[f_idx]]
            dr_idx = dd.dims.index('drift_rate')
            drate  = dd.drift_rate[idx[dr_idx]]
            data = blank_hit(data, metadata, f0, drate, padding=8)
            hit_id += 1
        else:
            break
    return hits

hits = find_hits_recursive(frame.data)

Then plotting:

data = np.copy(np.expand_dims(frame.data, axis=1))
data = normalize(data.astype('float32'), return_space='cpu')
dd, _md = dedoppler(data, metadata, max_dd=1.0, return_space='cpu')
plt.imshow(dd.data.squeeze(), aspect='auto')
hits = find_hits_recursive(frame.data)
h = np.array(hits)
plt.scatter(h[:, 2], h[:, 0], color='red', marker='x')

image

Main advantage is that it can find hits hiding behind very bright signals. Would likely need a lot of optimizations if used in practice, but could be good to get rid of high S/N hits as a first-pass

telegraphic commented 2 years ago

image

telegraphic commented 2 years ago

Now implemented 💥