Power spectrum of downsampled image

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

• What is your criteria for saying that one downsampled spectrum "follows" or not follows the original closely? It seems to me that the 3 are perfectly valid representations of the red one. They all have some mean squared error. It seems on the lower left image, the donwsampled spectrum has a tiny bit of a bias (it looks to be on average more below than above), but that shouldn't be caused by downsampling. – bone Sep 24 '15 at 13:18
• This is just a by-eye observation. As you say, the bottom-left is below the red spectrum - significantly below it might add, rather than a little bit. The bottom two pictures should show qualitatively the same thing, but they clearly do not. I agree that I wouldn't know why or how downsampling causes this, but to my mind the downsampling is the only difference between the two bottom panels and so the error HAS to be there. – user1991 Sep 24 '15 at 13:50