Electroencephalogram (EEG) signals may get easily contaminated by muscle artifacts, which may lead to wrong interpretation in the brain-computer interface (BCI) system as well as in various medical diagnoses. The main objective of this paper is to remove muscle artifacts without distorting the information contained in the EEG. A novel multi-stage EEG denoising method is proposed for the first time where wavelet packet decomposition (WPD) is combined with modified non-local means (NLM) algorithm. At first, the artifacted EEG is identified through a pre-trained classifier. Next, the identified EEG signal is decomposed into wavelet coefficients through WPD. Muscle artifacts are eliminated from the wavelet coefficients by estimating the clean wavelet coefficients through a modified NLM algorithm instead of thresholding them. Finally, the artifact-free EEG is reconstructed from corrected wavelet coefficients through inverse WPD. To optimize the filter parameters of the NLM algorithm, two met-heuristic algorithms, grey wolf optimization (GWO) and particle swarm optimization (PSO), are used in this paper for the first time. The proposed system is first validated on simulated EEG data and then tested on real EEG data. The proposed approach achieved average mutual information (MI) as 2.9684 (0.7045) on real EEG data. The result reveals that the proposed system outperforms the recently developed denoising techniques with higher average MI, which indicates that the proposed approach is better in terms of quality of reconstruction and is fully automatic and can be implemented in online applications.