PowerCutStats/analysis.py at main · QTNM/PowerCutStats · GitHub

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
from __init__ import *


def get_signal_peaks(signal, cut):
    """_summary_

    Args:
        signal (times ): Time series of voltage data, the signal
        cut (_type_): Value upon which to cut FFT spectrum peaks

    Returns:
        _type_: _description_
    """

    signal_FFT = np.fft.fft(signal, norm='forward')

    abs_signal_FFT = np.abs(signal_FFT)

    peaks, _ = abs_signal_FFT[abs_signal_FFT > cut], np.where(abs_signal_FFT > cut)

    return peaks #, np.where(abs_signal_FFT > cut)


def get_dists(peaks, y, tau_1f, Nsamp):
    """_summary_

    Args:
        peaks (_type_): _description_

    Returns:
        _type_: _description_
    """

    if peaks == []:

        cdf = (1-np.exp(-(y**2)/tau_1f)) ** (Nsamp)
        pdf = np.gradient(cdf, y[1]-y[0])

        return cdf, pdf

    rice_cdf = np.ones(20001)
    for peak in peaks:
        rice_cdf *= scipy.stats.rice.cdf(y, b=abs(peak)/np.sqrt(tau_1f/2), loc=0, scale=np.sqrt(tau_1f/2))

    signal_cdf = rice_cdf * (1-np.exp(-(y**2)/tau_1f)) ** (Nsamp-peaks.size)

    signal_pdf = np.gradient(signal_cdf, y[1]-y[0])

    return signal_cdf, signal_pdf

def calculate_noise_dists(vals, tau, Nsamp):
    ray_cdf = (1-np.exp(-(vals**2)/tau)) ** Nsamp
    ray_pdf = np.gradient(ray_cdf, vals[1]-vals[0])

    return ray_cdf, ray_pdf