Images from CIFAR-10 dataset

%reset -f
import time

start = time.time()

from IPython.display import display, HTML

display(HTML('''
<script id="MathJax-script"
        async
        src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-svg.js">
</script>
'''))

import warnings
import matplotlib.pyplot as plt
from pathlib import Path
import plotly.graph_objects as go
from torch.utils.data import DataLoader
from torchvision import datasets, transforms
import torch
import numpy as np
import plotly
from plotly.subplots import make_subplots
import os

from dyned import (
    sigma_delta_quantisation,
    generate_step_signal,
    identify_peaks_and_troughs,
    create_spike_plot,
    create_rle_spike_plot,
    create_lzw_spike_plot,
    create_huffman_spike_plot,
    create_delta_spike_plot,
    create_bwt_spike_plot,
    create_lz77_spike_plot,
    create_dynedc_v4_spike_plot,
    compression_summary,
    print_compression_table,
)

plotly.offline.init_notebook_mode(connected=True)

import plotly.io as pio
if os.environ.get("theme", "light") == "dark":
    pio.templates.default = "plotly_dark"
else:
    pio.templates.default = "plotly"

HERE = Path.cwd()

if torch.cuda.is_available():
    device = torch.device("cuda")
    print("Using GPU")
else:
    device = torch.device("cpu")
    print("Using CPU")

def get_mean_std(loader):
    # Var[X] = E[X^2] - E[X]^2
    channels_sum, channels_squared_sum, num_batches = 0, 0, 0

    for data, _ in loader:
        channels_sum += torch.mean(data, dim=[0, 2, 3])
        channels_squared_sum += torch.mean(data ** 2, dim=[0, 2, 3])
        num_batches += 1

    _mean = channels_sum / num_batches
    _std = (channels_squared_sum / num_batches - _mean ** 2) ** 0.5

    return _mean, _std

train_transform = transforms.Compose([transforms.ToTensor()])
train_dataset = datasets.CIFAR10(f"{HERE}/assets", train=True, download=True, transform=train_transform)
train_dataloader = DataLoader(train_dataset, batch_size=64, shuffle=True, num_workers=0)
train_mean, train_std = get_mean_std(train_dataloader)
print(f"Train Mean: {train_mean}, Train Standard Deviation: {train_std}")

test_transform = transforms.Compose([transforms.ToTensor()])
test_dataset = datasets.CIFAR10(f"{HERE}/assets", train=False, download=True, transform=test_transform)
test_dataloader = DataLoader(test_dataset, batch_size=64, shuffle=False, num_workers=0)
test_mean, test_std = get_mean_std(test_dataloader)
print(f"Test Mean: {test_mean}, Test Standard Deviation: {test_std}")


transform_train = transforms.Compose([transforms.ToTensor(), transforms.Normalize(train_mean.tolist(), train_std.tolist())])
train_dataset = datasets.CIFAR10(f"{HERE}/assets", train=True, download=True, transform=transform_train)

transform_test = transforms.Compose([transforms.ToTensor(), transforms.Normalize(test_mean.tolist(), test_std.tolist())])
test_dataset = datasets.CIFAR10(f"{HERE}/assets", train=False, download=True, transform=transform_test)

dataset = test_dataset

# # Find the mean and standard deviation of the ImageNet dataset
# transform = transforms.Compose([transforms.ToTensor()])
# dataset = datasets.ImageNet(f"{HERE}/assets", train=True, download=True, transform=transform, split="balanced")
# dataloader = DataLoader(dataset, batch_size=64, shuffle=False, num_workers=0)
# mean, std = get_mean_std(dataloader)
# print(f"Mean: {mean}, Standard Deviation: {std}")

# # Normalize the dataset
# transform = transforms.Compose([transforms.ToTensor(), transforms.Normalize(mean, std)])
# dataset = datasets.ImageNet(f"{HERE}/assets", train=True, download=True, transform=transform, split="balanced")
# dataloader = DataLoader(dataset, batch_size=64, shuffle=True, num_workers=4)

Using GPU

/home/akshay/code/personal/head-tracker/.venv/lib/python3.13/site-packages/torchvision/datasets/cifar.py:83: VisibleDeprecationWarning: dtype(): align should be passed as Python or NumPy boolean but got `align=0`. Did you mean to pass a tuple to create a subarray type? (Deprecated NumPy 2.4)
  entry = pickle.load(f, encoding="latin1")

Train Mean: tensor([0.4913, 0.4821, 0.4465]), Train Standard Deviation: tensor([0.2470, 0.2435, 0.2616])

Test Mean: tensor([0.4942, 0.4851, 0.4504]), Test Standard Deviation: tensor([0.2467, 0.2430, 0.2616])

img, label = dataset[3]
class_name = dataset.classes[label]

# 1. Cleaner PyTorch-native unnormalize function
def unnormalize(tensor, mean, std):
    """Reverse normalization strictly in PyTorch before returning."""
    # Reshape 1D mean/std tensors to (3, 1, 1) for broadcasting
    m = mean.view(3, 1, 1)
    s = std.view(3, 1, 1)

    # Unnormalize and clip to [0, 1]
    unnorm_tensor = torch.clamp(tensor * s + m, 0, 1)
    return unnorm_tensor

# 2. Unnormalize using the TRAIN mean and std
img_unnorm = unnormalize(img, train_mean, train_std)

# 3. Calculate grayscale on the UNNORMALIZED tensor
img_gray = 0.2989 * img_unnorm[0] + 0.5870 * img_unnorm[1] + 0.1140 * img_unnorm[2]

# 4. Prepare for Matplotlib (convert to NumPy and permute channels to H, W, C)
img_rgb = img_unnorm.permute(1, 2, 0).numpy()

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(3, 1.5), dpi=150) # Slightly larger for better visibility
fig.suptitle(f"Label: {label}, Class: {class_name}", fontsize=8)

axes[0].imshow(img_rgb, interpolation='bilinear')
axes[0].set_title("RGB", fontsize=7)
axes[0].axis('off')

axes[1].imshow(img_gray.numpy(), cmap='gray', interpolation='bilinear')
axes[1].set_title("Grayscale", fontsize=7)
axes[1].axis('off')

plt.tight_layout()
plt.show()

# Flatten the image and compute FFT
flatten_img = torch.flatten(img_gray)
fft_result = torch.fft.fft2(img_gray)
fft_shifted = torch.fft.fftshift(fft_result)
fft_shifted_flatten = torch.flatten(fft_shifted)
# magnitude_spectrum = 20 * torch.log10(torch.abs(fft_shifted_flatten) + 1e-6)
magnitude_spectrum = torch.abs(fft_shifted_flatten)
phase_spectrum = torch.angle(fft_shifted_flatten)

magnitude_spectrum = torch.log1p(magnitude_spectrum)
magnitude_spectrum_np = magnitude_spectrum.numpy()

# Calculate frequency bins
N = len(magnitude_spectrum_np)
sample_rate = 1 # Unknown sample rate, so we assume it is 1 Hz
freqs = np.fft.fftshift(np.fft.fftfreq(N, d=1/sample_rate))

# Create Plotly subplots
fig = make_subplots(rows=2, cols=1, subplot_titles=("Flattened Image Data", "1D FFT Magnitude Spectrum"), vertical_spacing=0.3)

# Add flattened image data
fig.add_trace(
    go.Scatter(x=np.arange(len(flatten_img)), y=flatten_img.numpy(), mode='lines', name='Flattened Image Data'),
    row=1, col=1
)
fig.update_xaxes(title_text="Index", row=1, col=1)
fig.update_yaxes(title_text="Pixel Value", row=1, col=1)

# Add FFT magnitude spectrum
fig.add_trace(
    go.Scatter(x=freqs, y=magnitude_spectrum, mode='lines', name='Magnitude Spectrum', line=dict(color='blue')),
    row=2, col=1
)
fig.update_xaxes(title_text="Frequency (Hz)", row=2, col=1)
fig.update_yaxes(title_text="Magnitude (dB)", row=2, col=1)

fig.update_xaxes(title_text="Frequency (Hz)", row=3, col=1)
fig.update_yaxes(title_text="Magnitude (dB) Peaks", row=3, col=1)

fig.update_layout(height=500, width=1400, title_text=f"Analysis of CIFAR-10 Class: {class_name}", showlegend=False)
fig.write_html(f"{HERE}/{class_name}_{label}_EMNIST_FFT.html")
fig.show()

# Quantise the magnitude spectrum using DyNED encoder
quantised_magnitude_spectrum, errors = sigma_delta_quantisation(magnitude_spectrum_np, levels=20)
step_signal = generate_step_signal(quantised_magnitude_spectrum)
peaks, troughs = identify_peaks_and_troughs(quantised_magnitude_spectrum)

fig = make_subplots(rows=5, cols=1, vertical_spacing=0.09, subplot_titles=(
    "Magnitude Spectrum", "Quantised Signals", "Spikes", "Peaks and Troughs", "Quantisation Error"))

fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum_np, name="Magnitude Spectrum"), row=1, col=1)

fig.add_trace(go.Scatter(x=freqs, y=quantised_magnitude_spectrum, mode='lines', name='Quantised Magnitude Spectrum'), row=2, col=1)

fig.add_trace(create_spike_plot(freqs, step_signal), row=3, col=1)

fig.add_trace(go.Scatter(x=freqs[peaks], y=quantised_magnitude_spectrum[peaks], mode='markers', marker=dict(color='green', size=3), name='Peaks'), row=4, col=1)
fig.add_trace(go.Scatter(x=freqs[troughs], y=quantised_magnitude_spectrum[troughs], mode='markers', marker=dict(color='red', size=3), name='Troughs'), row=4, col=1)

fig.add_trace(go.Scatter(x=np.arange(len(errors)), y=errors, mode='lines', name='Quantisation Errors', line=dict(color='red')), row=5, col=1)

fig.update_layout(height=1080, width=1400, title='Quantised FFT Magnitude Spectrum', xaxis_title='Frequency (Hz)', yaxis_title='Magnitude (dB)')
fig.show()

Number of spikes: 732

Reconstruction of the Image

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from scipy.interpolate import RBFInterpolator, CubicSpline, Rbf, KroghInterpolator

def split_data(x, y, num_segments):
    """Split the data into approximately equal segments."""
    # Determine the size of each segment
    segment_size = len(x) // num_segments
    segments = []

    for i in range(num_segments):
        start = i * segment_size
        # Ensure the last segment grabs all remaining data
        if i == num_segments - 1:
            end = len(x)
        else:
            end = start + segment_size
        segments.append((x[start:end], y[start:end]))

    return segments

def interpolate_segments(segments):
    """Interpolate each segment and return interpolators."""
    interpolators = []
    for x_seg, y_seg in segments:
        interpolator = KroghInterpolator(x_seg, y_seg)
        interpolators.append(interpolator)
    return interpolators

def merge_interpolations(interpolators, x_full):
    """Evaluate interpolations on the full x range, merging them smoothly."""
    y_full = np.zeros_like(x_full)
    segment_length = len(x_full) // len(interpolators)

    for i, interpolator in enumerate(interpolators):
        if i == len(interpolators) - 1:
            # Last segment goes up to the end
            x_range = x_full[i * segment_length:]
        else:
            x_range = x_full[i * segment_length: (i + 1) * segment_length]
        y_full[i * segment_length: i * segment_length + len(x_range)] = interpolator(x_range)

    return y_full

# Assuming quantised_magnitude_spectrum is already defined
quantised_magnitude_spectrum_reshaped = quantised_magnitude_spectrum.reshape(-1, 1)
freqs_reshaped = freqs.reshape(-1, 1)  # Ensure this is a 2D array as required by RBFInterpolator

freqs_interpolated = np.linspace(freqs.min(), freqs.max(), len(quantised_magnitude_spectrum))  # Smooth interpolation across the frequency range

# RBF interpolation
start_time = time.time()
rbfi = RBFInterpolator(freqs_reshaped, quantised_magnitude_spectrum_reshaped, kernel="linear")
interpolated_values = rbfi(freqs_interpolated[:, None]).flatten()
rbf_time = time.time() - start_time

# Cubic spline interpolation
start_time = time.time()
cube_spline = CubicSpline(freqs, quantised_magnitude_spectrum, bc_type='clamped')
interpolated_values_cubic = cube_spline(freqs_interpolated)
cubic_time = time.time() - start_time

# Radial basis function interpolation
start_time = time.time()
rbf = Rbf(freqs, quantised_magnitude_spectrum, function='cubic')
interpolated_values_rbf = rbf(freqs_interpolated)
rbf_classic_time = time.time() - start_time

# Krogh interpolation
start_time = time.time()
segments = split_data(freqs, quantised_magnitude_spectrum, num_segments=30)
interpolators = interpolate_segments(segments)
interpolated_values_krogh = merge_interpolations(interpolators, freqs)
krogh_time = time.time() - start_time

# Plot using Plotly
fig = make_subplots(rows=4, cols=1, subplot_titles=("RBF Interpolation", "Cubic Interpolation", "RBF Interpolation (Classic)", "Krogh Interpolation"), vertical_spacing=0.09)
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum, mode='lines', name='Quantized Magnitude Spectrum'), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values, mode='lines',line = dict(dash='dot'), name='Interpolated Quantized Magnitude Spectrum'), row=1, col=1)

fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum, mode='lines', name='Quantized Magnitude Spectrum'), row=2, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_cubic, mode='lines', line = dict(dash='dot'), name='Interpolated Quantized Magnitude Spectrum'), row=2, col=1)

fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum, mode='lines', name='Quantized Magnitude Spectrum'), row=3, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_rbf, mode='lines', line = dict(dash='dot'), name='Interpolated Quantized Magnitude Spectrum'), row=3, col=1)

fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum, mode='lines', name='Quantized Magnitude Spectrum'), row=4, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_krogh, mode='lines', line = dict(dash='dot'), name='Interpolated Quantized Magnitude Spectrum'), row=4, col=1)
fig.update_layout(title='Interpolated Quantized FFT Magnitude Spectrum', xaxis_title='Frequency (Hz)', yaxis_title='Magnitude (dB)', height=1080, width=1400, showlegend=False)
fig.show()

original_values_at_interpolated_freqs = np.interp(freqs_interpolated, freqs, magnitude_spectrum)

# Compute interpolation method
# RBF Interpolation Error Calculation
mse_rbfi = mean_squared_error(original_values_at_interpolated_freqs, interpolated_values)
rmse_rbfi = np.sqrt(mse_rbfi)
mae_rbfi = mean_absolute_error(original_values_at_interpolated_freqs, interpolated_values)
r2_rbfi = r2_score(original_values_at_interpolated_freqs, interpolated_values)

# Cubic Spline Interpolation Error Calculation
mse_cubic = mean_squared_error(original_values_at_interpolated_freqs, interpolated_values_cubic)
rmse_cubic = np.sqrt(mse_cubic)
mae_cubic = mean_absolute_error(original_values_at_interpolated_freqs, interpolated_values_cubic)
r2_cubic = r2_score(original_values_at_interpolated_freqs, interpolated_values_cubic)

# Radial Basis Function Interpolation Error Calculation
mse_rbf = mean_squared_error(original_values_at_interpolated_freqs, interpolated_values_rbf)
rmse_rbf = np.sqrt(mse_rbf)
mae_rbf = mean_absolute_error(original_values_at_interpolated_freqs, interpolated_values_rbf)
r2_rbf = r2_score(original_values_at_interpolated_freqs, interpolated_values_rbf)

# Krogh Interpolation Error Calculation
mse_krogh = mean_squared_error(original_values_at_interpolated_freqs, interpolated_values_krogh)
rmse_krogh = np.sqrt(mse_krogh)
mae_krogh = mean_absolute_error(original_values_at_interpolated_freqs, interpolated_values_krogh)
r2_krogh = r2_score(original_values_at_interpolated_freqs, interpolated_values_krogh)

# Print Error Metrics
print(f"RBFInterpolator - RMSE: {rmse_rbfi}, MAE: {mae_rbfi}, R2: {r2_rbfi}, Time: {rbf_time:.4f}s")
print(f"Cubic Spline - RMSE: {rmse_cubic}, MAE: {mae_cubic}, R2: {r2_cubic}, Time: {cubic_time:.4f}s")
print(f"Radial Basis Function - RMSE: {rmse_rbf}, MAE: {mae_rbf}, R2: {r2_rbf}, Time: {rbf_classic_time:.4f}s")
print(f"Krogh Interpolator - RMSE: {rmse_krogh}, MAE: {mae_krogh}, R2: {r2_krogh}, Time: {krogh_time:.4f}s")

fig = make_subplots(rows=1, cols=1)
fig.add_trace(go.Scatter(x=freqs, y=magnitude_spectrum, mode='lines', name='Original Data Points'), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values, mode='lines', name='RBF Interpolation', line=dict(dash='dot')), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_cubic, name='Cubic Spline', line=dict(dash='dot')), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_rbf, name='RBF', line=dict(dash='dot')), row=1, col=1)
fig.add_trace(go.Scatter(x=freqs_interpolated, y=interpolated_values_krogh, name='Krogh', line=dict(dash='dot')), row=1, col=1)
fig.update_layout(title='Comparison of Interpolation Methods', xaxis_title='Frequency (Hz)', yaxis_title='Magnitude (dB)', height=500, width=1400)
fig.show()

/tmp/ipykernel_2067169/1621631898.py:25: UserWarning: 34 degrees provided, degrees higher than about thirty cause problems with numerical instability with 'KroghInterpolator'
  interpolator = KroghInterpolator(x_seg, y_seg)
/tmp/ipykernel_2067169/1621631898.py:25: UserWarning: 38 degrees provided, degrees higher than about thirty cause problems with numerical instability with 'KroghInterpolator'
  interpolator = KroghInterpolator(x_seg, y_seg)

RBFInterpolator - RMSE: 0.13884961178192629, MAE: 0.11323501172777736, R2: 0.9642931858134081, Time: 0.0230s
Cubic Spline - RMSE: 0.13884961178159766, MAE: 0.1132350117275947, R2: 0.964293185813577, Time: 0.0005s
Radial Basis Function - RMSE: 0.13884955534723653, MAE: 0.11323491069651936, R2: 0.9642932148390934, Time: 0.0864s
Krogh Interpolator - RMSE: 0.13888957497235452, MAE: 0.11335826740113117, R2: 0.9642726288446478, Time: 0.0395s

def reconstruct_image(interpolated_magnitude, original_phase, original_shape):
    # Combine interpolated magnitude with original phase
    complex_spectrum = interpolated_magnitude * torch.exp(1j * original_phase)

    # Reshape back to 2D
    complex_spectrum_2d = complex_spectrum.reshape(original_shape)

    # Perform inverse FFT shift and then inverse FFT
    ifft_result = torch.fft.ifft2(torch.fft.ifftshift(complex_spectrum_2d))

    # Take the real part
    reconstructed_image = torch.real(ifft_result)

    return reconstructed_image

interpolated_magnitude_2d = interpolated_values_cubic.reshape(magnitude_spectrum.shape)
reconstructed_img = reconstruct_image(torch.exp(torch.from_numpy(interpolated_magnitude_2d)) - 1, phase_spectrum, img_gray.shape)

img_np = img_gray.cpu().numpy()          # (32, 32) grayscale
magnitude_spectrum_np = magnitude_spectrum.cpu().numpy()
reconstructed_img_np = reconstructed_img.cpu().numpy()

fig, axes = plt.subplots(1, 2, figsize=(3, 1.5), dpi=150)
fig.suptitle("Image Reconstruction", fontsize=7)

axes[0].imshow(img_np, cmap='gray', interpolation='bilinear')
axes[0].set_title("Original Image", fontsize=5)
axes[0].axis('off')

axes[1].imshow(reconstructed_img_np, cmap='gray', interpolation='bilinear')
axes[1].set_title("Reconstructed Image", fontsize=5)
axes[1].axis('off')

plt.tight_layout()
plt.show()

Structural Similarity Index (SSI)

from skimage.metrics import structural_similarity as ssim

similarity_index = ssim(img_np, reconstructed_img_np, data_range=reconstructed_img_np.max() - reconstructed_img_np.min())
print(f"Structural Similarity Index: {similarity_index}")

Structural Similarity Index: 0.9841686132181577

Convolutions

import torch
from torch.nn import Conv2d
import numpy as np

def apply_convolution(_img, _kernel):
    if not isinstance(_img, torch.Tensor):
        _img = torch.tensor(_img, dtype=torch.float32)
    if _img.dim() == 2:
        _img = _img.unsqueeze(0).unsqueeze(0)
    elif _img.dim() == 3:
        _img = _img.unsqueeze(0)

    conv_layer = Conv2d(in_channels=1, out_channels=1, kernel_size=3, padding=1, bias=False)
    conv_layer.weight.data = _kernel.reshape(1, 1, 3, 3)

    with torch.no_grad():
        output = conv_layer(_img)

    return output.squeeze().numpy()

def visualize_multiple_convolutions(original, images: dict[str, np.ndarray]):
    n = len(images) + 1
    cols = 3
    rows = (n + cols - 1) // cols

    fig, axes = plt.subplots(rows, cols, figsize=(4, rows * 1.3), dpi=150)
    axes = axes.flatten()

    axes[0].imshow(original, cmap='gray', interpolation='bilinear')
    axes[0].set_title("Original Image", fontsize=5)
    axes[0].axis('off')

    for i, (name, _img) in enumerate(images.items()):
        axes[i + 1].imshow(_img, cmap='gray', interpolation='bilinear')
        axes[i + 1].set_title(name, fontsize=5)
        axes[i + 1].axis('off')

    # Hide unused subplots
    for j in range(n, len(axes)):
        axes[j].axis('off')

    fig.suptitle("Multiple Convolutions Visualization", fontsize=7)
    plt.tight_layout()
    plt.show()

# Define 5 different kernels
kernels = {
    "Laplacian": torch.tensor([[0, 1, 0], [1, -4, 1], [0, 1, 0]], dtype=torch.float32),  # Laplacian
    "Sharpen": torch.tensor([[-1, -1, -1], [-1, 8, -1], [-1, -1, -1]], dtype=torch.float32),  # Sharpen
    "Gaussian Blur": torch.tensor([[1, 2, 1], [2, 4, 2], [1, 2, 1]], dtype=torch.float32) / 16,  # Gaussian Blur
    "Sobel X": torch.tensor([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype=torch.float32),  # Sobel X
    "Sobel Y": torch.tensor([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype=torch.float32) # Sobel Y
}

convolved_images = {name: apply_convolution(img_np, kernel) for name, kernel in kernels.items()}
visualize_multiple_convolutions(img_np, convolved_images)

# Print statistics
print(f"Original Image - Min: {img.min():.4f}, Max: {img.max():.4f}, Mean: {img.mean():.4f}")
for i, (name, convolved) in enumerate(convolved_images.items()):
    print(f"Convolution using {name} - Min: {convolved.min():.4f}, Max: {convolved.max():.4f}, Mean: {convolved.mean():.4f}")

Original Image - Min: -1.4641, Max: 1.7859, Mean: 0.5171
Convolution using Laplacian - Min: -1.3483, Max: 0.6665, Mean: -0.0691
Convolution using Sharpen - Min: -1.7263, Max: 3.3967, Mean: 0.2054
Convolution using Gaussian Blur - Min: 0.1440, Max: 0.8453, Mean: 0.5771
Convolution using Sobel X - Min: -3.3525, Max: 3.3836, Mean: 0.0036
Convolution using Sobel Y - Min: -2.1869, Max: 3.2668, Mean: -0.0558

end = time.time()
print(f"Time taken: {end - start:.4f} seconds")

Time taken: 8.8956 seconds

Spike Compression

# Compression classes are imported from dyned.py
# (RunLengthEncoding, LZWCompressor, HuffmanCompressor, DeltaCompressor,
#  BWTransform, LZ77Compressor, DyNEDcCompressor)
# See imports in the first cell.

plt1, cr1, sa1 = create_rle_spike_plot(freqs, step_signal)
plt2, cr2, sa2 = create_lzw_spike_plot(freqs, step_signal)
plt3, cr3, sa3 = create_huffman_spike_plot(freqs, step_signal)
plt4, cr4, sa4 = create_delta_spike_plot(freqs, step_signal)
plt5, cr5, sa5 = create_bwt_spike_plot(freqs, step_signal)
plt6, cr6, sa6 = create_lz77_spike_plot(freqs, step_signal)
plt7, cr7, sa7 = create_dynedc_v4_spike_plot(freqs, step_signal)

fig = make_subplots(rows=8, cols=1, subplot_titles=(
    f"Original Spikes",
    f"Run-Length Encoding ({cr1}, {sa1})",
    f"LZW Compression ({cr2}, {sa2})",
    f"Huffman Compression ({cr3}, {sa3})",
    f"Delta Compression ({cr4}, {sa4})",
    f"BWT Compression ({cr5}, {sa5})",
    f"LZ77 Compression ({cr6}, {sa6})",
    f"$\\text{{DyNED}}_{{c}}\\text{{ Compression ({cr7}, {sa7})}}$"
    ), vertical_spacing=0.09)

fig.add_trace(create_spike_plot(freqs, step_signal), row=1, col=1)

fig.add_trace(plt1, row=2, col=1)
fig.add_trace(plt2, row=3, col=1)
fig.add_trace(plt3, row=4, col=1)
fig.add_trace(plt4, row=5, col=1)
fig.add_trace(plt5, row=6, col=1)
fig.add_trace(plt6, row=7, col=1)
fig.add_trace(plt7, row=8, col=1)

fig.update_layout(height=1080, width=1400, title='Spike Compression Techniques - <i>r</i> image', xaxis_title='Frequency (Hz)', yaxis_title='Magnitude (dB)', showlegend=False)
fig.show()

# Compression summary table with bits-per-spike comparison
summary = compression_summary(step_signal)
print_compression_table(summary)

Number of spikes: 732
Run-Length Encoding:
    Original data size: 1024 bits
    Compressed data size: 2217 bits
    Compression ratio: 2.165
    Space saving: -116.5%
Number of spikes: 732
LZW Compression:
    Original data size: 1024 bits
    Compressed data size: 1456 bits
    Compression ratio: 1.422
    Space saving: -42.2%
Number of spikes: 732
Huffman Compression:
    Original data size: 1024 bits
    Compressed data size: 871 bits
    Compression ratio: 0.851
    Space saving: 14.9%
Number of spikes: 732
Delta Compression:
    Original data size: 1024 bits
    Compressed data size: 1680 bits
    Compression ratio: 1.641
    Space saving: -64.1%
Number of spikes: 732
Burrows-Wheeler Transform:
    Original data size: 1024 bits
    Compressed data size: 1031 bits
    Compression ratio: 1.007
    Space saving: -0.7%
Number of spikes: 732
LZ77:
    Original data size: 1024 bits
    Compressed data size: 2409 bits
    Compression ratio: 2.353
    Space saving: -135.3%
Number of spikes: 732
    Mode selected: huffman
DyNEDc Compression:
    Original data size: 1024 bits
    Compressed data size: 871 bits
    Compression ratio: 0.851
    Space saving: 14.9%
Number of spikes: 732

Spike Train: 1024 bits, 732 spikes
Method         Compressed    Ratio   Bits/Spike   Saving
--------------------------------------------------------
RLE                  2217   2.1650         3.03  -116.5%
LZW                  1456   1.4219         1.99   -42.2%
Huffman               871   0.8506         1.19    14.9%
Delta                1680   1.6406         2.30   -64.1%
BWT                  1031   1.0068         1.41    -0.7%
LZ77                 2409   2.3525         3.29  -135.3%
DyNEDc                951   0.9287         1.30     7.1%
DyNEDc V2            1110   1.0840         1.52    -8.4%
DyNEDc V3             871   0.8506         1.19    14.9%
DyNEDc V4             871   0.8506         1.19    14.9%