Positron emission tomography (PET) images suffer from a noticeable amount of statistical noise. In order to reduce the noise, a post smoothing process is usually adopted in the conventional PET systems. However, its performance is limited because the process is mostly based on a Gaussian random noise which is quite distinct from the noise of PET images. It has been reported that noise variance of each voxel is proportional to the square of the mean value in a PET image reconstructed by expectation-maximization (EM). In addition, we also observe that the variance varies with the spatial sensitivity distribution in a PET system. Based on those properties, we determine a unique formula representing a relationship between the mean and variance for a given PET system. Meanwhile, a block matching 3D (or 4D) algorithm is known as the state of the art in Gaussian noise reduction. To effectively apply it for noise reduction of PET images, we first perform a noise characteristic conversion from the PET image noise to Gaussian random noise using a pre-determined relationship. We then apply a block matching 4D (BM4D) algorithm and reconvert the result. Using the Monte Carlo simulation, we demonstrate that proposed algorithm can effectively reduce the noise in the whole image region while minimizing the image resolution degradation.