An active area of neuroimaging research involves examining functional relationships between spatially remote brain regions. When determining whether two brain regions exhibit significant correlation due to true functional connectivity, one must account for the background spatial correlation inherent in neuroimaging data. We define background correlation as spatiotemporal correlation in the data caused by factors other than neurophysiologically based functional associations such as scanner induced correlations and image preprocessing. We develop a 4D spatiotemporal wavelet packet resampling method which generates surrogate data that preserves only the average background spatial correlation within an axial slice, across axial slices, and through each voxel time series, while excluding the specific correlations due to true functional relationships. We also extend an amplitude adjustment algorithm which adjusts our surrogate data to closely match the amplitude distribution of the original data. Our method improves upon existing wavelet-based methods and extends them to 4D. We apply our resampling technique to determine significant functional connectivity from resting state and motor task fMRI datasets.