QFC - 用于电生理学的时序数据量化傅里叶压缩

概述

随着细胞外电生理学数据的规模不断增加，开发压缩多通道时间序列数据的有效方法至关重要。虽然无损方法是理想的选择，可以完美保留原始信号，但这些方法的压缩比通常只有2-4倍。我们需要的是10-30倍的压缩比，这促使我们考虑有损方法。

在这里，我们实现了一种简单的有损压缩方法，该方法受到离散余弦变换（DCT）和JPEG压缩图像的量化步骤的启发。该方法包括以下步骤

• 计算时间域中时序数据的离散傅里叶变换（DFT）。
• 量化傅里叶系数以实现目标熵（熵决定了理论上可实现的压缩比）。这是通过乘以归一化因子然后四舍五入到最接近的整数来完成的。
• 使用zlib或zstd（可以使用其他方法）压缩减少熵的量化傅里叶系数。

解压缩

• 解压缩量化傅里叶系数。
• 除以归一化因子。
• 计算逆离散傅里叶变换（IDFT）以获得重建的时间序列数据。

这种方法特别适用于经过带通滤波的数据，因为抑制的傅里叶系数会产生量化信号的特别低熵。

有关各种有损和无损压缩方案的比较，请参阅大规模电生理学数据的压缩策略，Buccino等。。

安装

pip install qfc

示例用法

# See examples/example1.py

from matplotlib import pyplot as plt
import numpy as np
from qfc import qfc_estimate_quant_scale_factor
from qfc.codecs import QFCCodec


def main():
    sampling_frequency = 30000
    duration = 2
    num_channels = 10
    num_samples = int(sampling_frequency * duration)
    y = np.random.randn(num_samples, num_channels) * 50
    y = lowpass_filter(y, sampling_frequency, 6000)
    y = np.ascontiguousarray(y)  # compressor requires C-order arrays
    y = y.astype(np.int16)
    target_residual_stdev = 5

    ############################################################
    quant_scale_factor = qfc_estimate_quant_scale_factor(
        y,
        target_residual_stdev=target_residual_stdev
    )
    codec = QFCCodec(
        quant_scale_factor=quant_scale_factor,
        dtype="int16",
        segment_length=10000,
        compression_method="zstd",
        zstd_level=3
    )
    compressed_bytes = codec.encode(y)
    y_reconstructed = codec.decode(compressed_bytes)
    ############################################################

    y_resid = y - y_reconstructed
    original_size = y.nbytes
    compressed_size = len(compressed_bytes)
    compression_ratio = original_size / compressed_size
    print(f"Original size: {original_size} bytes")
    print(f"Compressed size: {compressed_size} bytes")
    print(f"Actual compression ratio: {compression_ratio}")
    print(f'Target residual std. dev.: {target_residual_stdev:.2f}')
    print(f'Actual Std. dev. of residual: {np.std(y_resid):.2f}')

    xgrid = np.arange(y.shape[0]) / sampling_frequency
    ch = 3  # select a channel to plot
    n = 1000  # number of samples to plot
    plt.figure()
    plt.plot(xgrid[:n], y[:n, ch], label="Original")
    plt.plot(xgrid[:n], y_reconstructed[:n, ch], label="Decompressed")
    plt.plot(xgrid[:n], y_resid[:n, ch], label="Residual")
    plt.xlabel("Time")
    plt.title(f'QFC compression ratio: {compression_ratio:.2f}')
    plt.legend()
    plt.show()


def lowpass_filter(input_array, sampling_frequency, cutoff_frequency):
    F = np.fft.fft(input_array, axis=0)
    N = input_array.shape[0]
    freqs = np.fft.fftfreq(N, d=1 / sampling_frequency)
    sigma = cutoff_frequency / 3
    window = np.exp(-np.square(freqs) / (2 * sigma**2))
    F_filtered = F * window[:, None]
    filtered_array = np.fft.ifft(F_filtered, axis=0)
    return np.real(filtered_array)


if __name__ == "__main__":
    main()

Zarr示例

见examples/zarr_example.py

基准测试

我已经在真实和合成数据上建立了一些初步的系统基准测试。请参见./benchmarks和./benchmarks/results。

如所见

• 在傅里叶域量化（QFC）比在时域量化（称为QTC）对于真实数据或带通滤波数据要好得多。
• 与未滤波的原始数据相比，带通滤波数据的压缩比要好得多。
• 对于方法的无损部分，zstd比zlib好，在压缩比、压缩速度和解压速度这三个因素上都更好。
• 显然，压缩比将严重依赖于目标残差标准差。

许可证

此代码根据Apache许可证版本2.0提供。

作者

Jeremy Magland，计算数学中心，Flatiron研究所