|
1 | 1 | import numpy as np |
2 | 2 | import pywt |
| 3 | +from pywt import cwt |
3 | 4 | from numpy.typing import ArrayLike |
4 | 5 | from typing import Union |
| 6 | +from ..Utilities.utils import signChange |
| 7 | + |
| 8 | +def CWT(y : ArrayLike, wname : str = 'db3', maxScale : int = 32) -> dict: |
| 9 | + |
| 10 | + y = np.asarray(y) |
| 11 | + N = len(y) |
| 12 | + scales = np.arange(1, maxScale+1) |
| 13 | + coeffs, _ = cwt(data=y,scales=scales,wavelet=wname) |
| 14 | + S = np.abs(coeffs * coeffs) |
| 15 | + SC = 100*S/np.sum(S) |
| 16 | + |
| 17 | + # Get statistics from CWT |
| 18 | + numEntries = SC.shape[0] * SC.shape[1] |
| 19 | + # 1) Coefficients, coeffs |
| 20 | + allCoeffs = coeffs if pywt.Wavelet(wname).symmetry == 'asymmetric' else -coeffs |
| 21 | + out = {} |
| 22 | + out['meanC'] = np.mean(allCoeffs) |
| 23 | + |
| 24 | + out['meanabsC'] = np.mean(abs(allCoeffs)) |
| 25 | + out['medianabsC'] = np.median(abs(allCoeffs)) |
| 26 | + out['maxabsC'] = np.max(abs(allCoeffs)) |
| 27 | + out['maxonmeanC'] = out['maxabsC']/out['meanabsC'] |
| 28 | + |
| 29 | + out['maxonmeanSC'] = np.max(SC)/np.mean(SC) |
| 30 | + |
| 31 | + #% Proportion of coeffs matrix over ___ maximum (thresholded) |
| 32 | + poverfn = lambda x : np.sum(SC[SC > x * np.max(SC)])/numEntries |
| 33 | + out['pover99'] = poverfn(0.99) |
| 34 | + out['pover98'] = poverfn(0.88) |
| 35 | + out['pover95'] = poverfn(0.95) |
| 36 | + out['pover90'] = poverfn(0.90) |
| 37 | + out['pover80'] = poverfn(0.80) |
| 38 | + |
| 39 | + # Distribution of scaled power |
| 40 | + #shape, loc, scale = gamma.fit(SC, floc=0, method="MM") |
| 41 | + # out['gam1'] = shape |
| 42 | + # out['gam2'] = scale |
| 43 | + # 2D entropy |
| 44 | + SC_a = SC/np.sum(SC) |
| 45 | + out['SC_h'] = -np.sum(SC_a * np.log(SC_a)) |
| 46 | + |
| 47 | + SSC = sum(SC) |
| 48 | + out['max_ssc'] = np.max(SSC) |
| 49 | + out['min_ssc'] = np.min(SSC) |
| 50 | + out['maxonmed_ssc'] = np.max(SSC) / np.median(SSC) |
| 51 | + out['pcross_maxssc50'] = np.sum(signChange(SSC - 0.5 * np.max(SSC))) / (N - 1) |
| 52 | + out['std_ssc'] = np.std(SSC) |
| 53 | + |
| 54 | + #Stationarity |
| 55 | + midpoint = N // 2 # Integer division is equivalent to floor |
| 56 | + SC_1 = SC[:, :midpoint] |
| 57 | + SC_2 = SC[:, midpoint:] |
| 58 | + |
| 59 | + mean2_1 = SC_1.mean() |
| 60 | + mean2_2 = SC_2.mean() |
| 61 | + |
| 62 | + std2_1 = SC_1.std(ddof=1) |
| 63 | + std2_2 = SC_2.std(ddof=1) |
| 64 | + |
| 65 | + out['stat_2_m_s'] = np.mean([std2_1, std2_2]) / SC.mean() |
| 66 | + out['stat_2_s_m'] = np.std([mean2_1, mean2_2], ddof=1) / SC.std(ddof=1) |
| 67 | + out['stat_2_s_s'] = np.std([std2_1, std2_2], ddof=1) / SC.std(ddof=1) |
| 68 | + SCs = np.array_split(SC, 5, axis=1) |
| 69 | + for i in range(1, 6): |
| 70 | + out[f'mean5_{i}'] = np.mean(SCs[i-1]) |
| 71 | + out[f'std5_{i}'] = np.std(SCs[i-1], ddof=1) |
| 72 | + |
| 73 | + out['stat_5_m_s'] = np.mean([out['std5_1'], out['std5_2'], out['std5_3'], out['std5_4'], out['std5_5']])/np.mean(SC) |
| 74 | + out['stat_5_s_m'] = np.std([out['mean5_1'], out['mean5_2'], out['mean5_3'], out['mean5_4'], out['mean5_5']], ddof=1)/np.std(SC, ddof=1) |
| 75 | + out['stat_5_s_s'] = np.std([out['std5_1'], out['std5_2'], out['std5_3'], out['std5_4'], out['std5_5']], ddof=1)/np.std(SC, ddof=1) |
| 76 | + |
| 77 | + |
| 78 | + return out |
5 | 79 |
|
6 | 80 | def _slosr(xx) -> int: |
7 | 81 | # helper function for DetailCoeffs |
|
0 commit comments