update test1.py
1 changed file+89−11
test1.pymodified+89−11View file
@@ -1,19 +1,25 @@
11 # %%
22 import numpy as np
3-from zia._filters import bandpass_filter
3+from zia._filters import bandpass_filter, highpass_filter
44 from zia._data_loaders import load_real_000876, load_real_000409, load_real_001290
55 from zia._compress_ints_lossless import compress_ints_lossless
66 from zia._analysis import linear_fit, compute_entropy_per_sample, estimate_noise_level
7+import matplotlib.pyplot as plt
78
89 # %%
910 N = 500_000
1011
11-# X = load_real_001290(num_samples=N, num_channels=1, start_channel=0)
12-X = load_real_000409(num_samples=N, num_channels=1, start_channel=101)
13-# X = load_real_000876(num_samples=N, num_channels=1, start_channel=45)
12+channel_number = 101
13+X = load_real_000409(num_samples=N, num_channels=1, start_channel=channel_number).flatten()
14+
15+# X = load_real_001290(num_samples=N, num_channels=1, start_channel=0).flatten()
16+# X = load_real_000876(num_samples=N, num_channels=1, start_channel=45).flatten()
1417
1518 X = X.astype(np.int16)
16-X = X.flatten()
19+
20+# %%
21+plt.figure(figsize=(12, 4))
22+plt.plot(X[:2400])
1723 # %%
1824 def print_ideal_compression_ratio(X):
1925 ee = compute_entropy_per_sample(X)
@@ -63,11 +69,16 @@ print_ideal_compression_ratio(X_mr)
6369 v = 0.25 # step size for quantization
6470 lowcut = 300
6571 highcut = 6000
66-X2 = bandpass_filter(X - np.median(X), sampling_frequency=30000, lowcut=lowcut, highcut=highcut)
67-noise_level = estimate_noise_level(X2, sampling_frequency=30000)
68-X2b = X2 / noise_level / v
72+X_filt = bandpass_filter(X - np.median(X), sampling_frequency=30000, lowcut=lowcut, highcut=highcut)
73+noise_level = estimate_noise_level(X_filt, sampling_frequency=30000)
74+X_filt_normalized = X_filt / noise_level
75+X2b = X_filt_normalized / v
6976 X2 = np.round(X2b).astype(np.int16)
7077
78+# %%
79+plt.figure(figsize=(12, 4))
80+plt.plot(X[:2400])
81+
7182 # %%
7283 print('FILTERED (and quantized)')
7384 print_ideal_compression_ratio(X2)
@@ -84,15 +95,82 @@ print_ideal_compression_ratio(residuals2)
8495 # %%
8596 import matplotlib.pyplot as plt
8697 plt.figure(figsize=(12, 4))
87-plt.plot(X[:800])
98+plt.plot(X[:2400])
8899 plt.title('RAW')
89100
90101 plt.figure(figsize=(12, 4))
91-plt.plot(X2[:800])
102+plt.plot(X2[:2400])
92103 plt.title('FILTERED')
93104
94105 plt.figure(figsize=(12, 4))
95-plt.plot(residuals2[:800])
106+plt.plot(residuals2[:2400])
96107 plt.title('FILTERED MARCOVIAN')
97108
98109 # %%
110+def sliding_max(x, delta):
111+ y = np.zeros_like(x)
112+ for i in range(len(x)):
113+ y[i] = np.max(x[max(0, i - delta):min(len(x), i + delta + 1)])
114+ return y
115+
116+def smoothed(x, delta):
117+ y = np.zeros_like(x)
118+ for i in range(len(x)):
119+ y[i] = np.mean(x[max(0, i - delta):min(len(x), i + delta + 1)])
120+ return y
121+
122+cc = [3, 6]
123+Y = sliding_max(np.abs(X_filt_normalized), 50)
124+Y = smoothed(Y, 20)
125+Y = np.minimum(1, np.maximum(0, (Y - cc[0]) / (cc[1] - cc[0])))
126+# Y = highpass_filter(Y, sampling_frequency=30000, lowcut=3)
127+Y_scaled = X_filt_normalized * Y
128+X3b = Y_scaled / v
129+X3 = np.round(X3b).astype(np.int16)
130+
131+plt.figure(figsize=(12, 4))
132+plt.plot(X2[:2400], color='lightgray')
133+# plt.plot(Y[4000:5000])
134+plt.plot(X3b[:2400])
135+# %%
136+print('FILTERED (and quantized) with suppression')
137+print_ideal_compression_ratio(X3)
138+print('')
139+print('FILTERED DELTA ENCODING with suppression')
140+print_ideal_compression_ratio(np.diff(X3))
141+print('')
142+print('FILTERED MARCOVIAN with suppression')
143+residuals3 = get_marcovian_prediction_residual(X3, 20)
144+print_ideal_compression_ratio(residuals3)
145+print('ACTUAL FILTERED MARCOVIAN with suppression')
146+print_actual_compression_ratios(residuals3)
147+# %%
148+def get_run_lengths(x):
149+ runs = []
150+ i = 0
151+ current_nonzero_run_length = 0
152+ while i < len(x):
153+ if np.all(x[i:i+10] == 0):
154+ runs.append(current_nonzero_run_length)
155+ current_nonzero_run_length = 0
156+ j = i
157+ while j < len(x) and x[j] == 0:
158+ j += 1
159+ runs.append(j - i)
160+ i = j
161+ else:
162+ current_nonzero_run_length += 1
163+ i += 1
164+ if np.max(runs) < 256:
165+ return np.array(runs, dtype=np.uint8)
166+ if np.max(runs) < 2 ** 16:
167+ return np.array(runs, dtype=np.uint16)
168+ return np.array(runs, dtype=np.uint32)
169+
170+AA = residuals3[residuals3 != 0]
171+run_lengths = get_run_lengths(residuals3)
172+print(run_lengths, run_lengths.itemsize)
173+ee = compute_entropy_per_sample(AA)
174+theoretical_size = (len(AA) * ee / 8) + run_lengths.nbytes
175+theoretical_compression_ratio = len(residuals3) * X.itemsize / theoretical_size
176+print(f'Theoretical compression ratio: {theoretical_compression_ratio:.2f}')