Freqbot¶

import numpy as np
import scipy.stats, scipy.signal
from PIL import ImageStat
from matplotlib.pylab import find

def is_greyscale(im):
    """
    Decides if an imae is greyscale.
    """
    stat = ImageStat.Stat(im)
    if sum(stat.sum[:3])/3 == stat.sum[0]:
        return True
    return False

def parabolic(f, x):
    """
    Interpolation.
    """
    xv = 1/2. * (f[x-1] - f[x+1]) / (f[x-1] - 2 * f[x] + f[x+1]) + x
    yv = f[x] - 1/4. * (f[x-1] - f[x+1]) * (xv - x)
    return (xv, yv)

def freq_from_crossings(sig, fs):
    """
    Frequency from zero crossings.
    """
    indices = find((sig[1:] >= 0) & (sig[:-1] < 0))
    crossings = [i - sig[i] / (sig[i+1] - sig[i]) for i in indices]
    return fs / np.mean(np.diff(crossings))

def freq_from_autocorr(sig, fs):
    """
    Frequency from autocorrelation.
    """
    sig = sig + 128
    corr = scipy.signal.fftconvolve(sig, sig[::-1], mode='full')
    corr = corr[len(corr)/2:]
    d = np.diff(corr)
    start = find(d > 0)[0]
    peak = np.argmax(corr[start:]) + start
    px, py = parabolic(corr, peak)
    return fs / px

def get_spectrum(signal, fs):
    """
    Computes the spectrum. Returns frequencies, coefficients,
    f_min and f_max.
    """
    windowed = signal * scipy.signal.blackmanharris(len(signal))
    a = abs(np.fft.rfft(windowed))
    f = np.fft.rfftfreq(len(signal), 1/fs)

    db = 20 * np.log10(a)
    sig = db - np.amax(db) + 20
    indices = find((sig[1:] >= 0) & (sig[:-1] < 0))
    crossings = [z - sig[z] / (sig[z+1] - sig[z]) for z in indices]
    mi, ma = np.amin(crossings), np.amax(crossings)
    x = np.arange(0, len(f))  # for back-interpolation
    f_min = np.interp(mi, x, f)
    f_max = np.interp(ma, x, f)
    return f, a, f_min, f_max

def freq_from_fft(signal, fs):
    """
    Peak frequency from FFT.
    """
    f, a, f_min, f_max = get_spectrum(signal, fs)
    i = np.argmax(a)
    true_i = parabolic(np.log(a), i)[0]
    return fs * true_i / len(signal)

def get_snr(i):
    """
    Bad algorithm
    """
    i += 128
    return np.nanmean(i) / np.nanstd(i)

def get_phase(i):
    """
    Computes phase using method from Purves's tutorial in TLE.
    """
    e = scipy.signal.hilbert(i)

    # Get the biggest 25 sample indices and sort them by amplitude
    biggest = np.argpartition(e.real, -25)[-25:]
    s = np.vstack((biggest, e.real[biggest])).T
    sort = s[s[:,1].argsort()][::-1]

    # Prune the list down to the biggest for realz
    biggest_pruned = [sort[:,0][0]]
    for ix in sort[:,0][1:]:
        add = True
        for got in biggest_pruned:
            if abs(ix - got) < 5:  # made-up number
                add = False
        if add:
            biggest_pruned.append(ix)
            if len(biggest_pruned) == 5: break

    # Get the interpolated phase values
    results = []
    for ix in biggest_pruned:
        true_i = parabolic(np.log(abs(e)+0.01), ix)[0]
        x = np.arange(0, len(e))
        rad = np.interp(true_i, x, np.angle(e))
        results.append(np.degrees(rad))

    return np.nanmean(results)

def get_trace_indices(y, ntraces, spacing):
    """
    Returns the indices of the traces to analyse.
    """
    if spacing == 'random':
        x = 0.05 + 0.9*np.random.random(ntraces)  # avoids edges
        ti = np.sort(x * y)
    else:
        n = ntraces + 1
        ti = np.arange(1./n,1.,1./n) * y
    return np.round(ti).astype(int)

def analyse(i, t_min, t_max, trace_indices, func):
    """
    Does everything.
    """
    fs = i.shape[0] / (t_max - t_min)

    spec, freq, phase, snr = [], [], [], []
    mis, mas = [], []
    for ti in trace_indices:
        trace = i[:, ti]
        try:
            f = func(trace, fs)
            freq.append(f)

            p = get_phase(trace)
            phase.append(p)
            
            snr.append(get_snr(trace))

            frq, amp, fmi, fma = get_spectrum(trace, fs)
            spec.append(amp)
            mis.append(fmi)
            mas.append(fma)

        except Exception:
            continue

    return spec, freq, phase, snr, mis, mas

sig = np.array([1, 2, 3, 2, 1, 0, -1, -2, -1, 1, 3, 1, -1, -3, 1, 2])
i1 = ((sig[1:] >= 0) & (sig[:-1] < 0)).nonzero()
i2 = find((sig[1:] >= 0) & (sig[:-1] < 0))

# url = "https://math.berkeley.edu/~sethian/2006/Applications/Seismic/time_mig_img.jpg"
url = "http://www.spectrumgeo.com/wp-content/uploads/BPC-Bahamas-Data-Example.jpg"

method = 'xing'
avg = 'mean'
region = None
# region = '199,141,1608,1087'
ntraces = 10
trace_spacing = 'random' or 'regular'
spectrum = 'none'
bins = 9
t_min = 0
t_max = 3.4

width, height = im.size

height = 50

dts = [0.0005, 0.001, 0.002, 0.004, 0.008]

for dt in sorted(dts, reverse=True):
    target = (t_max - t_min) / dt
    if target > height: break