1
$\begingroup$

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?

enter image description here

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
|improve this question
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – bone Sep 24 '15 at 13:18
  • $\begingroup$ 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. $\endgroup$ – user1991 Sep 24 '15 at 13:50
1
$\begingroup$

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:

enter image description here

|improve this answer
$\endgroup$
  • $\begingroup$ Nice of you to share your solution. $\endgroup$ – David Aug 19 '18 at 17:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.