-
Notifications
You must be signed in to change notification settings - Fork 0
/
testing.py
154 lines (116 loc) · 3.97 KB
/
testing.py
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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
"""CVM
from CVM import cvm
from CVM_python import cvm as cvm_python
# Check results
print("Input z:")
z = int(input())
print("Input win_len:")
win_len = int(input())
print("Cython", cvm(z, win_len))
print("Python", cvm_python(z, win_len))
"""
"""SNR
from SNR import snr
from SNR_python import snr as snr_python
# Check results
x = [0,1,2,3]
y = [4,5,6,7]
print("Cython", snr(x, y))
print("Python", snr_python(x, y))
"""
"""MSE
from MSE import mse
from MSE_python import mse as mse_python
# Check results
x = [0,1,2,3,8,6,4,1,2,9,7,3,5,5,6,2,1,9,0,8]
y = [4,5,6,7,1,8,7,2,9,8,4,2,6,8,8,4,1,6,4,7]
print("Cython", mse(x, y))
print("Python", mse_python(x, y))
"""
"""CDFCALC
from CDFCALC import cdfcalc
from CDFCALC_python import cdfcalc as cdfcalc_python
# Check results
x = [0,1,2,3,8,6,4,1,2,9,7]
disn = 7
ind = 4
print("Cython", list(cdfcalc(x, disn, ind)))
print("Python", cdfcalc_python(x, disn, ind))
"""
"""ecdf
from ecdf import ecdf
from ecdf_python import ecdf as ecdf_python
# Check results
array = [101, 118, 121, 103, 142, 111, 119, 122, 128, 112, 117,157]
y, x = ecdf(array)
y_p, x_p = ecdf_python(array)
print("Cython\n", f"\nx: {list(x)}\ny: {list(y)}")
print("\nPython\n", f"\nx: {x_p}\ny: {y_p}")
"""
"""VMD
import numpy as np
import xarray as xr
from Prop_VMD_CVM import Prop_VMD_CVM
from Prop_VMD_CVM_python import Prop_VMD_CVM as Prop_VMD_CVM_python
# Input parameter
win_len = 32 # Window length
NIMF = 10 # IMFs to consider
Np = 36 # No of consecutive iterations that must detect for the detection to hold
N_mon = 10
sigL = 12
# Select Pfa using a decaying function e^(-k+1)
opt_Pfa = np.exp(-np.arange(0, NIMF, 1))
# Load signal
a = xr.load_dataset('Data/2015_Granada_Noisy.nc').beta_mean.values[0]
a1 = xr.load_dataset('Data/2015_Granada_Denoised.nc').beta_mean.values[0]
# Generating the noisy signal for given input SNR
iSNR = 20
SNRii = iSNR
sigma = 0.1 + (0.3 / (SNRii + 1))
x = np.arange(-1/2, 1/2, 1/(win_len-1))
g = np.exp(- x**2 / (2 * sigma**2))
g /= sum(g)
snri = iSNR
f = a
# Proposed approach
imf, rec, y = Prop_VMD_CVM(a, f, win_len, NIMF, opt_Pfa, Np)
imf_p, rec_p, y_p = Prop_VMD_CVM_python(a, f, win_len, NIMF, opt_Pfa, Np)
# Check results are the same
print(f"Are imfs the same: {sum(sum(imf == imf_p))/(imf_p.shape[0] * imf_p.shape[1]) == True}")
print(f"Are recs the same: {sum(sum(rec == rec_p))/(rec_p.shape[0] * rec_p.shape[1]) == True}")
print(f"Are ys the same: {sum(sum(y == y_p))/(y_p.shape[0] * y_p.shape[1]) == True}")
"""
import numpy as np
import xarray as xr
from vmd_cvm_cython.Prop_VMD_CVM import Prop_VMD_CVM
from vmd_cvm_python.Prop_VMD_CVM_python import Prop_VMD_CVM as Prop_VMD_CVM_python
import matplotlib.pyplot as plt
# Input parameter
win_len = 32 # Window length
NIMF = 10 # IMFs to consider
Np = 36 # No of consecutive iterations that must detect for the detection to hold
# Select Pfa using a decaying function e^(-k+1)
opt_Pfa = np.exp(-np.arange(0, NIMF, 1))
# Load signal
a = xr.load_dataset('Data/2015_Granada_Noisy.nc').beta_mean.values[0].astype(np.float64)
range_ = xr.load_dataset('Data/2015_Granada_Noisy.nc')['range']
f = a
# Proposed approach
imf, imf_hat, omega = Prop_VMD_CVM(a, f, win_len, NIMF, opt_Pfa, Np)
imf_p, imf_hat_p, omega_p = Prop_VMD_CVM_python(a, f, win_len, NIMF, opt_Pfa, Np)
# Plot results
# fig, ax = plt.subplots(1, 2, figsize=(10, 6), sharey=True)
# ax[0].plot(sum(imf), range_, color='cornflowerblue', linewidth=0.7, label='Cython')
# ax[0].set_xlabel('signal [arb. units]')
# ax[0].set_ylabel('range [m]')
# ax[0].set_title('Cython reconstruction')
# ax[0].grid()
# ax[1].plot(sum(imf_p), range_, color='forestgreen', linewidth=0.7, label='Python')
# ax[1].set_xlabel('signal [arb. units]')
# ax[1].set_title('Python reconstruction')
# ax[1].grid()
# plt.show()
# Check results are the same
print(f"Are imfs the same: {imf.shape == imf_p.shape}")
print(f"Are recs the same: {imf_hat.shape == imf_hat_p.shape}")
print(f"Are ys the same: {omega.shape == omega_p.shape}")