aboutsummaryrefslogtreecommitdiffstats
path: root/src/two_tone_lib.py
blob: df7d53f468e4fdc06279adc8abdb2464c25b8495 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
import numpy as np
import matplotlib.pyplot as plt

def gen_two_tone(path = "./input.dat", predist = None, par = None, debug = False):
    period1 = 3.875
    period2 = 4
    t_both = 124
    assert(t_both / period1 % 1 == 0)
    assert(t_both / period2 % 1 == 0)

    t = np.arange(0,t_both)
    sin1 = np.sin(t * 2 * np.pi * 1./period1)
    sin2 = np.sin(t * 2 * np.pi * 1./period2)
    sig = sin1 + sin2

    if predist is None:
        res = sig
    else:
        res = predist(sig, par)

    res = res / np.max(res)

    res = res.astype(np.complex64)
    res.tofile(path)

    a_load = np.fromfile(path, dtype=np.complex64)
    assert(np.isclose(a_load, res).all()), "Inconsistent stored file"

    if debug == True:
        plt.plot(np.abs(np.concatenate((a_load, a_load))))
        plt.savefig(path + ".png")
        plt.clf()

        plt.plot(np.abs(np.fft.fftshift(np.fft.fft(np.concatenate((a_load, a_load))))), 'ro')
        plt.savefig(path + "_fft.png")
        plt.clf()

    return path

def predist_poly(sig, coefs = []):
    res = sig
    for idx, coef in enumerate(coefs):
        res += sig * np.abs(sig)**(idx+1) * coef #+1 because first correction term is squared
    return res

def analyse_power_spec(spec, debug = False, debug_path="", suffix=""):
    peak_1 = None
    peak_2 = None
    spec_start = 4096
    spec_end = 8192
    first_peak = spec_start + 2048
    second_peak = spec_start + 2114
    delta_freq = 66
    peak_other = []
    if debug: plt.plot(spec[spec_start:spec_end])
    for x in [c * delta_freq + delta_freq//2 for c in range(spec_start//delta_freq)]:
        start = spec_start + x 
        end = spec_start + x + delta_freq
        peak = spec[start:end].max()
        if debug: plt.plot((start-spec_start,end-spec_start), (peak, peak))
        if start < first_peak and end > first_peak:
            peak_1 = peak
            if debug: plt.plot((start-spec_start,end-spec_start), (peak+1, peak+1))
        elif start < second_peak and end > second_peak:
            peak_2 = peak
            if debug: plt.plot((start-spec_start,end-spec_start), (peak+1, peak+1))
        else:
            peak_other.append(peak)
    mean_signal = (peak_1 + peak_2) / 2
    mean_others = np.mean(peak_other)
    score = mean_signal - mean_others
    if debug:
        plt.savefig(debug_path + "/" + str(score) + suffix + ".png")
        plt.clf()
    return score