Power spectrum of downsampled image

Question

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

Answer 1

I have solved the problem.

Since I was taking the sum of 4 pixels in downsampling, what I was doing in effect was convolving with a sum of delta functions (that is, one delta function for each pixel). In Fourier space this translates to multiplying with a cosine, which explains the increasingly strong (relative) discrepancy between input and observed power spectrum.

To mitigate this, Wikipedia helped me out. If we first implement a lowpass filter in Fourier space (filtering out high frequencies, to avoid aliasing later on), and then reduce the sampling rate (to 1 in 4 pixels), we get the desired result, and indeed: the observed power spectrum looks nice up until the cutoff: