Functional magnetic resonance imaging (fMRI) is a major non-invasive method to detect brain activities with high spatial resolution. However, fMRI suffers from wide ranges of systematic and physiological noise leading to low functional signal to noise ratio [1]. To minimize artifacts, researchers have mostly focused on magnitude-based filtering methods and have largely neglected phase data due to its noisy nature [2]. However, fMRI is natively a complex-valued data with also useful phase information [3, 4]. Here, we propose a novel spatial weighted averaging method which uses the phase information along with magnitude to create a reference signal and utilize it to update weights iteratively. We evaluate the method on Averaged-BOSS (A-BOSS) fMRI data-set [5, 6] which has valuable phase information due to the inherent significant changes of phase within the transition band of Steady State Free Precession (SSFP) profile. The results indicate that the approach can suppress noise effectively and preserve the boundaries of active regions.Here, we propose a model-based anisotropic approach using phase information in addition to magnitude data. First, initial weights are calculated based on the Correlation Coefficient (CC) between the magnitude Time-Course (TC) of any voxel and its adjacent voxels. After phase unwrapping and masking data with phase quality map based on phase derivative variance algorithm [7], Fast Fourier Transform (FFT) is applied to complex TCs. Next, the initial weights and low-pass filtering are applied to the complex data. To estimate the reference, the real part of the signal in the time domain is chosen. In the following step, the similarity index is calculated based on squared CC between the magnitude TCs and the reference signal in the order to update weights. This procedure is repeated iteratively to reach convergence. The block diagram is shown in Fig. 1. To evaluate the method, we performed an experiment on a 3T Tim Trio Siemens scanner using TrueFISP. The scans were acquired at high spatial resolution (1.2*1.2*4 mm3). The method was applied on a healthy volunteer in a right hand finger tapping experiment in which 3 blocks of rest and act were used (rest/act=15s/15s). Other imaging parameters were: TR=30ms and FA=4°.No pre-processing steps such as motion correction, registration or cluster thresholding were performed for final results. This method is applied on an experimental data acquired by A-BOSS imaging technique. The magnitude and phase activation maps are shown respectively in Figs 2.a and 2b. Magnitude and phase data was processed based on correlation analysis of each TC with stimulus pattern and a p-value threshold of p=0.05 was chosen to detect active voxels. For the evaluation, we compared proposed method with Gaussian smoothing and magnitude-based (typical) anisotropic approach. Figs 3.a, 3.b and 3.c show the activity maps of Gaussian smoothing (kernel size equal to double size of each voxel), typical anisotropic smoothing, and proposed method (with 5 iterations) respectively. As depicted in activation maps, false positive (FP) voxels are eliminated successfully in our method compared to others and the active regions related to motor regions continuously grow and reasonably fit into the gray matter.The main issue in spatial filters is to eliminate the FP errors and preserve the boundaries of the active regions [8, 9]. Conventional anisotropic averaging is used based on spatial correlation of active neighboring voxels mostly in magnitude data. However, according to Fig 2.b, the phase activity map, also contains valuable spatial information. Therefore, we have exploited phase information along with magnitude data to create more accurate reference signal to improve the weights updating procedure. As shown in Fig. 3.a, the correlation map without any smoothing suffers from FPs. In Fig 3.b, Gaussian kernel has removed FPs but the active regions diffuse to neighboring regions. Although the typical anisotropic smoothing boundaries are preserved, this method suffers from FPs and non-continuous activated regions. In the proposed method, smoothing kernel minimizes the FPs, preserves the boundaries and also continuity of active regions. Furthermore, as shown in Fig 3.c, activation map fitting nicely into the gray matter at high spatial resolution. To evaluate the method numerically, highly activated voxels surviving the threshold of 67% of maximum correlation coefficient are selected in each map and the mean averages of them are shown in Table 1. As depicted in this table, the mean CC of Region of Interest (ROI) in our method is 81.31% which indicates the better performance compared to other methods. The proposed filter is also applicable for resting-state fMRI to improve the signal quality as well as task fMRI due to its independency to the stimulus pattern.