Skip to content
Snippets Groups Projects
Select Git revision
  • b73bc61cb38f4aab936a91f780f7036a402acf5e
  • main default protected
2 results

artefact_4_5_PSD.py

Blame
  • artefact_4_5_PSD.py 2.04 KiB
    import numpy as np
    import math
    from scipy.fft import fft, ifft
    
    def indices(a, func):
        return [i for (i, val) in enumerate(a) if func(val)]
    
    def nextpow2(x):
        """ Function for finding the next power of 2 """
        return 1 if x == 0 else math.ceil(math.log2(x))
    
    def powerspectrum(data, Fs):
        T = 1 / Fs
        L = np.size(data, 0)
        NFFT = 2 ** nextpow2(L)
        Y = fft(data, NFFT) / L
        f = Fs / 2 * np.linspace(0, 1, int(NFFT / 2))
        power = 2 * np.abs(Y[0:int(NFFT / 2)])
        return f, power
    
    """ 4).4-5 Check PSDs + PSD of gamma frequency """
    def PSD(ICs_projections, FS, distance, thGamma):
        gammaPSDT = []
        diff = []
        marked_ICs = np.zeros(64)
        for i in range(np.size(ICs_projections, 2)):
            if (marked_ICs[i] <= 3):
                actual_projection = np.transpose(ICs_projections[:, :, i])
                pre_diff = []
                pre_gamma = []
                for k in range(np.size(actual_projection, 0)):
                    f, power = powerspectrum(actual_projection[k, :], FS)
                    actual_idealDistro1 = f[1:]
                    actual_idealDistro2 = np.divide(1, f[1:])
    
                    normedActual_idealDistro2 = np.divide(actual_idealDistro2, np.max(actual_idealDistro2))
                    normedPsCheck = np.divide(power[1:], np.max(power[1:]))
    
                    actual_diff = normedPsCheck - normedActual_idealDistro2
    
                    actual_specDistT = np.mean(math.sqrt(np.dot(actual_diff, actual_diff.T)))
                    pre_diff.append(actual_specDistT)
                    actual_gammaPSDT = np.mean(power[f > 30])
                    pre_gamma.append(actual_gammaPSDT)
    
                diff.append(np.mean(np.array(pre_diff)))
                gammaPSDT.append(np.mean(np.array(pre_gamma)))
            else:
                diff.append(-1)
                gammaPSDT.append(-1)
        gammaIC = np.array(gammaPSDT)
        diffIC = np.array(diff)
    
        for i in range(np.size(ICs_projections, 2)):
            if (diffIC[i] > distance):
                marked_ICs[i] = marked_ICs[i] + 1
            if (gammaIC[i] > thGamma):
                marked_ICs[i] = marked_ICs[i] + 1
    
        return marked_ICs