def measure_states(states, measurement_matrix, measurement_covariance):
    """
    Measure a list of states with a measurement matrix in the presence of
    measurement noise.

    Args:
        states (array): states to measure. Shape is NxSTATE_DIM.
        measurement_matrix (array): Each state in *states* is measured with this
            matrix. Should be MEAS_DIMxSTATE_DIM in shape.
        measurement_covariance (array): Measurement noise covariance. Should be
            MEAS_DIMxMEAS_DIM.

    Returns:
        (array): NxMEAS_DIM array of measurements.

    """
    # Sanitise input
    measurement_matrix = np.atleast_2d(measurement_matrix)
    measurement_covariance = np.atleast_2d(measurement_covariance)
    measurement_dim = measurement_matrix.shape[0]
    if measurement_covariance.shape != (measurement_dim, measurement_dim):
        raise ValueError(("Measurement matrix and covariance have inconsistent "
                          "shapes {} and {}").format(measurement_matrix.shape,
                                                     measurement_covariance.shape))
    states = np.atleast_2d(states)

    # Special case: no output
    if states.shape[0] == 0:
        return np.zeros((0, measurement_dim))

    # Measure states
    measurement_means = measurement_matrix.dot(states.T).T
    measurement_noises = np.random.multivariate_normal(
        mean=np.zeros(measurement_dim), cov=measurement_covariance,
        size=states.shape[0]
    )

    return measurement_means + measurement_noises
Measure a list of states with a measurement matrix in the presence of
    measurement noise.

    Args:
        states (array): states to measure. Shape is NxSTATE_DIM.
        measurement_matrix (array): Each state in *states* is measured with this
            matrix. Should be MEAS_DIMxSTATE_DIM in shape.
        measurement_covariance (array): Measurement noise covariance. Should be
            MEAS_DIMxMEAS_DIM.

    Returns:
        (array): NxMEAS_DIM array of measurements.