Spatial Transfer Function vs. nPCA

A demonstration that the spatial transfer function can be improved to ~80% at ~300 arcseconds by reducing the number of PCA components subtracted. However, image fidelity and noise backgrounds suffer, so these sorts of maps are probably less ideal for point source extraction.

Below:

0 PCA

3 PCA

10 PCA

15 PCA

Spatial Transfer Functions, revisit 4

Last report was a bit of a fiasco. There were problems all over the place I didn't understand. I still don't but I've fixed them. My best guess at this point is that a pass-by-reference led to an unacceptable modification of an image. That doesn't even make sense - there was no place it could have happened - but, there you have it.

So, going through the process step by step.

This is the effect of smoothing an image:

Note from the 3rd figure that 100% recovery isn't reached until ~700 arcseconds.


Next question: What is happening at large spatial scales in the flat-spectrum simulations?

No obvious problems there.

Hmm, no apparent problem here either, though one might ask why the two curves approach each other in sky05 (alpha=-0.5).

So it appears that the reason for the bump up at low frequencies (long wavelengths) must be because of edge effects. After much hassle, I've addressed that by cropping images.

Finally, the averaged results:

So we've got an Official Spatial Transfer Function.

However, of course, we must note that there is a dependence on the atmosphere amplitude to source amplitude ratio: it appears that large-scale structure is *easier* to recover when the atmosphere is at higher amplitude. This makes sense: it is easier to distinguish faint astrophysical signal from bright atmosphere in this case. The reason I didn't run simulations to test this more is that the S/N ratio on small scales becomes poor for the low astrophysical amplitudes.

STF measurement attempts, round 3?

Below are plots of the spatial transfer function for a variety of simulation parameters. Each line represents the median (i.e., mean of the center 2 of 4) spatial transfer function of 4 realizations of the listed parameters.

A few unfortunate details are apparent:

1. The STF is background dependent, though not strongly
   
2. The STF appears to never exceed 85%, and be more typically 75%. If this is not matched by point-source PSDs, an explanation for the flux discrepancy becomes possible
   
3. The 50% recovery point is much closer in than previously stated, at closer to 150" than 300".
   

Note for comparison that the recovery fraction for point sources is much more nearly 1, at least in the calibrator type maps. We need to wait for Jared's simulations with point sources in a full size map to confirm this, but it seems pretty likely that the STF is different for point sources and extended structures.

Examining bolocat

Just in case we were wondering, the V1 bolocat is completely inconsistent with a lognormal distribution, but is perfectly consistent (or... at least reasonably consistent...) with a power law distribution.

These plots show power-law fits (red) and lognormal fits (blue) to the data (black). It's pretty obvious that the lognormal is a bad fit, but in case you're unconvinced, the ks test for the "source flux" has a probability 1.6e-7, and it is the highest likelihood by 17 orders of magnitude out of the 4 flux types.

By contrast, the simulations are on average (though not uniformly) more consistent with lognormal than powerlaw distributions:

Even in those examples where the KS test is slightly more favorable for the powerlaw distribution, the lognormal is a pretty good fit, and the failure points for the two distributions are in about the same place. The smoothness of the lognormal distribution is required to reproduce the observed distribution.

Note that the first 4 plots are for the whole BGPS survey. What about an individual field? For obvious reasons, I choose l30 again.

This gets to be a little more interesting - apparently the "source flux" has a tendency to pick up the power-law distributed background structure, since it is consistent with a lognormal (but note that it is also consistent with a powerlaw! The ks test doesn't really say definitively which is better). Although the fits look bad at high flux, note that this is a log-log plot and therefore the difference in probability is rather small.

What does this all indicate? It's not entirely clear whether individual fields are genuinely more lognormally-distributed or whether the number statistics are just worse. However, even the source flux is consistent with a power-law, while many realizations of the simulations are not.

Therefore, we should perform the next logical test - add point sources drawn from a power-law distribution (and a log-normal distribution?) and see what bolocat retrieves. We can at least say now that the point source contribution cannot be ignored, since there is no power-law distribution that can reproduce the observed bolocat flux distribution.

Bolocat examination part 1

This figure shows the flux recovery of bolocat. The Y axis shows the Input flux (green) and the Filtered flux (blue/cyan/purple) with a few different large-scale filters. The red line is the 1-1 line. Obviously, bolocat doesn't recover everything that was there - we removed a lot of flux, so recovery fractions are in the few % range. However, bolocat does a much better job than the simple filter at recovering positive fluxes. Therefore, the filter function still isn't good enough.

Testing out analytic filter functions

I've attempted to model the spatial filter function as a gaussian (or PSF) plus an inverse gaussian. i.e., the high-spatial-frequency components are smoothed with the PSF, and the low spatial frequency components are convolved with a (1-gaussian) high-pass filter.

First, the mildly good news: With a 300" FWHM large-scale cutoff, the filter PSD reasonably resembles the iterative map PSD:

Luckily, the double-filter goes a very long way in explaining the scale-free flux loss. In the following diagram, I show the effect of the filter compared to the input map.

The filter only recovers about 75% of the flux at ANY wavenumber. The map does slightly worse at high frequencies, which I can't explain yet.

These show the recovery fraction of the iterative maps, a gaussian smoothing function with FWHM=33", and the mid-pass-filter. Map20 (no smooth) has a lot of additional "noise power" at high spatial frequencies; if it wasn't for the telescope filter function, we would apparently have pretty good high-frequency recovery. Hmph.

Note that map20 is higher than the filter at some intermediate frequencies, but quite a bit lower at higher frequencies. Also note the moderately poor agreement between the 'smoothed' and 'smoothed (theory)' lines.

Finally, look at the comparison between map20 and fiiltered. The agreement is not bad for positive points; filtered is apparently slightly higher but that can be adjusted. The problem: the filter forces some structures that are negative or zero to be positive. For example, look at the feature at 210,300 that is negative in Map20 but positive in Filtered. In the real (input) map, this feature is lower than its surroundings - it is legitimately negative.

Masked edge effects?

I raised the possibility that the ratty edges of a scan could affect the power spectrum measurements, leading to a scale-free power loss. This doesn't appear to happen... instead, the masking adds some small-scale power and apparently a small amount of large-scale power.

Spatial Transfer Functions

The majority of the past week has been dedicated to debugging; it looks like cross-scanned simulations finally work.

The plot below is a derivation of the spatial transfer function for a number of different intrinsic sky power-law power spectra.

Justifying the above plot is essential.

First, the very steep power-laws [-3 in the example below] show a recovery fraction >1. This is simply because their S/N was inadequate - the output power spectrum is nearly flat, but at a level higher than the sky.

Second, the most plausible power-laws [-1.5 in the example below] show pretty good recovery (90-95% over the relevant range):

There are some "white" power losses, particularly in the flatter power-spectra. My best guess is that this has something to do with the relative scales being offset from a mean of 1, but so far all tests to show that that is the cause have in fact shown no problems at all. What else could cause a scale-independent power loss?

Also, the flat power spectrum (and inverted) aren't quite flat because I impose a "galactic scale height" on them. Should I stop doing that?

Measuring preferred angles in the Herschel data on larger scales

This is something of a repeat of yesterday's exercise, but for Herschel data on larger scales, where filamentation is expected (at different angles?) on a different scale than the Galactic plane.

First example: L030 500 microns

There is evidently a preferred direction that is correlated down to the resolution of the map, though the preference is smallest at the smallest scales. This may simply be a statement that there are more sources along the galactic plane, though, since there really isn't any particularly obvious filamentation in the image.

In L59, on the other hand, there is at most a very weak preference except at the largest scales corresponding to the Galactic Plane. This is somewhat interesting because there IS obvious filamentation in the L59 image, but it does not have a preferred direction. Unfortunately, it is not obvious whether or where filamentation shows up in fourier space if it does not have a preferred direction. There is no excess at any spatial scale that I can pick out.

Hunting for preferred directions

In our last group meeting, we discussed simulations using filamentary structures. I've been trying to determine how best to generate random/artificial filamentary structures in images. The first step in that direction is coming up with a way to measure asymmetries and therefore preferred directions within a real map. In order to do this, I've had to develop a number of new tools in agpy for azimuthally binned radial profiles and radially-averaged azimuthal profiles (radialprofile.py).

Examining real maps is necessary because a simple sinusoidal dependence of the power introduces filamentation along the same directions at ALL scales, which is not obviously (or in some cases, obviously not) the correct solution. Filaments are observed on large scales, but sometimes there can be 'kinks' in opposite directions on small scales.

So, the map I picked to examine was one with the most flagrantly obvious filamentary structure in it: the Motte DR 21 MAMBO map. It was also a choice of convenience because I already had the data on my laptop....

The preferred direction is quite obvious in this map: there is a long filament going up and down the map. Therefore, the DC component should be substantially higher in one direction than the other.

In the power-spectral-density image, it is quite clear that there is a preferred direction, though it is not obvious that the fourier transform is rotated 90 degrees from the image. My fourier intuition somewhat fails me here.... I realize that a broad, smooth profile in real space should be narrow and highly peaked in fourier space, but I don't fully understand why it effectively spreads out in the perpendicular direction, which I have confirmed that it does with a simple experiment.

I also don't know what the Shah function is, but it implies a periodic dip in the image at every 1/5th of the image, or every 50 pixels.

These are the power spectra averaged over different angles as labeled. -15 corresponds to -15 to +15, 15 corresponds to 15 to 45, etc. The highly peaked 75-105 power spectrum shows the large-scale filamentary profile. I think the difference in the DC component is actually an artifact of the azimuthal binning process: each pixel can only be assigned one angle, so the DC value isn't included in all of them... I'll need to find a workaround for that because it's quite deceptive.

The more interesting way to view the data - and perhaps to analyze maps - is to take radial averages in some range of spatial scales and plot the azimuthal dependence. There is a clear sinusoid at large scales. The legend shows "spatial frequency" in 1/pixel units. The distribution becomes more even with angle and even changes preferred direction at smaller scales (higher frequencies).

Next step is testing different approaches. I think an added, steeper-power-law component would probably be the best way to start.

Another suggestion, courtesy Bruce Elmegreen, is to attempt this sort of asymmetric power law sampling in 3 dimensions (with only 1 or 2 dimensions asymmetric) and then projecting down onto two dimensions.