def METHOD_NAME(self):
    # TODO: Most of this stuff is duplicated in MeanEstimator - move up the hierarchy?
    n = self.src.n
    L = self.src.L
    _2L = 2 * self.src.L
    kernel = np.zeros((_2L, _2L, _2L, _2L, _2L, _2L), dtype=self.dtype)
    sq_filters_f = np.square(evaluate_src_filters_on_grid(self.src))
    for i in trange(0, n, self.batch_size):
        _range = np.arange(i, min(n, i + self.batch_size))
        pts_rot = rotated_grids(L, self.src.rotations[_range, :, :])
        weights = sq_filters_f[:, :, _range]
        weights *= self.src.amplitudes[_range] ** 2
        if L % 2 == 0:
            weights[0, :, :] = 0
            weights[:, 0, :] = 0
        # TODO: This is where this differs from MeanEstimator
        pts_rot = np.moveaxis(pts_rot[::-1], 1, 0).reshape(-1, 3, L**2)
        weights = weights.T.reshape((-1, L**2))
        batch_n = weights.shape[0]
        factors = np.zeros((batch_n, _2L, _2L, _2L), dtype=self.dtype)
        for j in range(batch_n):
            factors[j] = anufft(weights[j], pts_rot[j], (_2L, _2L, _2L), real=True)
        factors = Volume(factors).to_vec()
        kernel += vecmat_to_volmat(factors.T @ factors) / (n * L**8)
    # Ensure symmetric kernel
    kernel[0, :, :, :, :, :] = 0
    kernel[:, 0, :, :, :, :] = 0
    kernel[:, :, 0, :, :, :] = 0
    kernel[:, :, :, 0, :, :] = 0
    kernel[:, :, :, :, 0, :] = 0
    kernel[:, :, :, :, :, 0] = 0
    logger.info("Computing non-centered Fourier Transform")
    kernel = fft.mdim_ifftshift(kernel, range(0, 6))
    kernel_f = fft.fftn(kernel)
    # Kernel is always symmetric in spatial domain and therefore real in Fourier
    kernel_f = np.real(kernel_f)
    return FourierKernel(kernel_f)