计算grid spacing

方式1：傅里叶变换 (Fourier Translation)

首先将rate map从空间域转换成频率，然后计算出最大的6个主频率（分别对应grid pattern的6个主峰，因为grid的主峰会均匀地分布在二维空间中），计算波矢长度，再转换成周期并取平均。

计算最大的6个主频率：

def analyze_grid_cell_spacing(ratemap, n=6, space_range=(-2, 2)):

    # 计算空间步长
    space_extent = space_range[1] - space_range[0]

    # 计算傅立叶变换
    fft_result = np.fft.fftshift(np.fft.fft2(ratemap))
    power_spectrum = np.abs(fft_result) ** 2  # 频谱
    
    # 找到最大的六个主频率（排除直流分量）
    center = np.array(power_spectrum.shape) // 2
    power_spectrum[center[0], center[1]] = 0  # 忽略直流分量
    peak_positions = np.argpartition(power_spectrum.flatten(), -n)[-n:] 
    # 这里的第一个n表示我们要找到最大的n个元素（主频率），第二个-n表示我们取分区后最后n个（即最大的n个）元素的索引
    peak_positions = np.array(np.unravel_index(peak_positions, power_spectrum.shape)).T

    # 计算每个主频率的波数和周期
    spacings = []
    for peak_pos in peak_positions:
        kx = 2 * np.pi * (peak_pos[1] - center[1]) / space_extent  # 波数kx（rad/m）
        ky = 2 * np.pi * (peak_pos[0] - center[0]) / space_extent  # 波数ky（rad/m）
        k = np.sqrt(kx**2 + ky**2)  # 总波数（rad/m）
        spacing = 2 * np.pi / k if k != 0 else np.inf  # 周期（m）
        spacings.append(spacing)

    # 取平均周期
    spacing = np.mean(spacings)

    return spacing

结果如下：

方式2：自相关图 (Autocorrelation)

实验类文章中的常用计算方法，参考以下文章method部分：

Sargolini F, Fyhn M, Hafting T, et al. Conjunctive representation of position, direction, and velocity in entorhinal cortex[J]. Science, 2006, 312(5774): 758-762.

实现代码如下（代码参考自Mikail）：

  def calculate_grid_cell_periodicity_and_orientation(
        rate_map: np.ndarray,
    ) -> Tuple[float, float, np.ndarray]:
        """
        Given a rate map with suspected grid cells, compute the

        :param rate_map:
        :return:
        """

        # Copied from Mikail
        # 1. 计算二维自相关矩阵：衡量信号与自身的相似性
        autocorr = correlate2d(rate_map, rate_map, mode="full")
        # 2. 寻找局部峰值
        tmp = autocorr
        maxes = peak_local_max(tmp)
        # 3. 计算峰值相对于中心的偏移
        maxes_from_center = maxes - np.array(tmp.shape) // 2
        # 4. 计算峰值的距离并排序
        distances = np.linalg.norm(maxes_from_center, axis=1)
        sorted_indices = np.argsort(distances)
        # 5. 选择最近的六个峰值
        nearest_six = maxes_from_center[sorted_indices[1:7]]
        # 6. 计算周期和方向
        period, period_err = np.average(np.linalg.norm(nearest_six, axis=1)), np.std(
            np.linalg.norm(nearest_six, axis=1)
        )
        orientations = np.arctan2(nearest_six[:, 1], nearest_six[:, 0])
        return period, period_err, orientations

最后将period * space_extent / resolustion，转换成真正物理空间的spacing。

文章来源于互联网:计算grid spacing