2import numpy as np
4from zia._data_loaders import load_real_000876, load_real_000409, load_real_001290
5from zia._compress_ints_lossless import compress_ints_lossless
6from zia._analysis import linear_fit, compute_entropy_per_sample, estimate_noise_level
a81e75aInitial commitJeremy Magland 7
8# %%
9N = 500_000
11# X = load_real_001290(num_samples=N, num_channels=1, start_channel=0)
12X = 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)
15X = X.astype(np.int16)
16X = X.flatten()
17# %%
19 ee = compute_entropy_per_sample(X)
20 print(f'Ideal compression ratio: {X.itemsize * 8 / ee:.2f} ({ee:.2f} bits per sample)')
22def print_actual_compression_ratios(X):
23 buf_zstd = compress_ints_lossless(X, method='zstd')
24 buf_zlib = compress_ints_lossless(X, method='zlib')
25 buf_lzma = compress_ints_lossless(X, method='lzma')
26 buf_ans = compress_ints_lossless(X, method='simple_ans')
27 print(f'Zstd compression ratio: {len(X) * X.itemsize / len(buf_zstd):.2f}')
28 print(f'Zlib compression ratio: {len(X) * X.itemsize / len(buf_zlib):.2f}')
29 print(f'Lzma compression ratio: {len(X) * X.itemsize / len(buf_lzma):.2f}')
30 print(f'simple_ans compression ratio: {len(X) * X.itemsize / len(buf_ans):.2f}')
32def get_marcovian_prediction_residual(X, M):
33 sequences = np.array([X[i:i+M] for i in range(len(X) - 2 * M + 1)])
34 predictors = sequences[:, :M - 1]
35 target = sequences[:, M - 1]
37 coeffs, predict = linear_fit(predictors, target)
38 predictions = predict(predictors)
39 predictions = np.round(predictions)
40 residuals = target - predictions
41 residuals = residuals.astype(np.int16)
42 return residuals
a81e75aInitial commitJeremy Magland 43
44# %%
46print_ideal_compression_ratio(X)
a81e75aInitial commitJeremy Magland 47
48# %%
50print_ideal_compression_ratio(np.diff(X))
a81e75aInitial commitJeremy Magland 51
52# %%
54print_actual_compression_ratios(np.diff(X))
55print_ideal_compression_ratio(np.diff(X))
a81e75aInitial commitJeremy Magland 56
57# %%
59print('RAW MARCOVIAN')
60print_ideal_compression_ratio(X_mr)
62# %%
63v = 0.25 # step size for quantization
65highcut = 6000
66X2 = bandpass_filter(X - np.median(X), sampling_frequency=30000, lowcut=lowcut, highcut=highcut)
67noise_level = estimate_noise_level(X2, sampling_frequency=30000)
71# %%
73print_ideal_compression_ratio(X2)
a81e75aInitial commitJeremy Magland 74
75# %%
77print_ideal_compression_ratio(np.diff(X2))
81print('FILTERED MARCOVIAN')
82print_ideal_compression_ratio(residuals2)
a81e75aInitial commitJeremy Magland 83
84# %%
85import matplotlib.pyplot as plt
87plt.plot(X[:800])
88plt.title('RAW')
90plt.figure(figsize=(12, 4))
91plt.plot(X2[:800])
92plt.title('FILTERED')
94plt.figure(figsize=(12, 4))
95plt.plot(residuals2[:800])
96plt.title('FILTERED MARCOVIAN')