I am not sure if this is the right stackexchange community for this question, but here goes.
I am generating a Gaussian Random Field (GRF) of a pre-defined power spectrum. In the image below, the input power spectrum is shown (in each panel) as the solid, smooth red line.
As a sanity check, I get the power spectrum of the resulting realised GRF. In the image below, I show this power spectrum for the original, high resolution random field realization in the upper right panel in blue. As we can see, it follows the input spectrum closely, as required.
Now, for an application I want to downsample the original image (by a factor of 2; that is, every 2x2 block of pixels becomes one pixel in the new image). If I then take the power spectrum of this downsampled (and thus lower-resolution) image, it does not follow the input power spectrum anymore; this is shown in the lower-left panel.
My first instinct was that this might be a resolution effect. Therefore, I generated another realization of the GRF, but this time at a lower resolution (indeed the same resolution as the downsampled image). If it was a resolution effect, this power spectrum should display the same discrepancy. It does not, however: in the lower right panel you can see that this lower resoltuion spectrum also follows the input spectrum closely.
Which leads me to believe that the discrepancy is caused solely by the downsampling. This is puzzling to me however, since downsampling seems so straightforward (perhaps deceptively so). Below, I also give the (python) code used to downsample the image.
So, rephrasing my question for clarity: Why is the power spectrum of the downsampled image not in qualitative agreement with the original power spectrum?
def downsample(self,im, fact=4): """ Downsample an image. Default with 2x2 blocks, but this can be adjusted if desired Needed because first we generate at high res, then apply lensing, then downsample to Planck resolution """ assert isinstance(fact, int), type(fact) sx, sy = im.shape X, Y = np.ogrid[0:sx, 0:sy] regions = sy/fact * (X/fact) + Y/fact res = ndimage.sum(im, labels=regions, index=np.arange(regions.max()+1)) res.shape = (sx/fact, sy/fact) self.downsampled = res