Version 0.6 begun
I started on (what I hope will be) v0.6 tonight. New things: -Sigma-rejection flagging in the spatial domain -Inverse-variance weighting
MAD flagger
; The Mad Flagger; Flag based on the median average deviation within a spatial pixelfunction mad_flagger,data,inds,flags,nsig=nsig t = systime(/sec) f0 = total(where(flags)) if n_e(nsig) eq 0 then nsig = 3 newflags = flags mx=max(inds) vec3=fltarr(mx+1) h=histogram(inds,reverse_indices=ri,OMIN=om) for j=0L,n_elements(h)-1 do begin if ri[j+1] gt ri[j] then begin v_inds = [ri[ri[j]:ri[j+1]-1]] if n_e(v_inds) gt 2 then begin vec = data[v_inds]; vecmad = mad(vec) ; the MAD is WAY too small! I ended up rejecting 8% of points! vecmad = stddev(vec) vecmed = median(vec,/even) madreject = where( (vec gt vecmed + nsig*vecmad) or (vec lt vecmed - nsig*vecmad) ) if (n_e(madreject) gt 0 and total(madreject)) gt 0 then begin reject_inds = v_inds[madreject] newflags[reject_inds] = 1 endif endif endif endfor print,"MAD flagger took ",strc(systime(/sec)-t)," seconds and flagged ",$ strc(round(total(where(newflags)) - f0)),' points' return,newflagsend
Median drizzling
http://groups.google.com/group/comp.lang.idl-pvwave/browse_thread/thread/762770933591238d Someone going by the name 'wox' introduced me to a spectacularly simple drizzling algorithm for medianing. It will be a little less efficient because it can't take advantage of some 'partial sum' tricks that the average can, but it will median over ALL data points, which is an advantage.
Massive upgrades, new strategy?
A few major things today:
#. I renamed 'ts_to_map' to 'drizzle' and made 'ts_to_map' a wrapper that will either median combine by scans (if passed the right keywords, i.e. tstomapmedian and scans_info). This will reject cosmic rays much more effectively than previously. #. My extensive test runs on the 'weird' fields completed. Conclusions:
- Deconvolution is DEFINITELY responsibly for the weird fuzzy effects
- seen esp. in l=33, l=2.
- Even a descending iterator is ineffective at alleviating this
- problem What I still don't understand is WHY the deconvolver is behaving badly for some regions. It is extracting peaks that are too high.
- Not deconvolving works pretty well, so I'm not really concerned. I'm
- resetting the defaults to no deconvolution once my next test is done....
- Started up a new test run to compare median stacking with weighted average stacking
5x1?
060609_o17 is weird looking. It appears to be TWO 3x1s with an overlap, but it's only 1 observation. Any idea what happened?
Weighting - next steps
Idea: 'edge taper' for scans. The scan start/end usually has higher noise, so it would be good to downweight those regions slightly. Probably only the first/last 5% of each scan should be linearly downweighted. I'm running a lot of tests on l001 with/without deconvolution, lots of iterations, different pca components. They'll be done sometime tomorrow.
Weighting and high-frequency noise
Image of PSDs (with no normalization) of the raw (blue), delined and exponential and polynomial subtracted (white), the noise timestream (yellow), and the data (cyan). The good: It looks like all of the powerline noise got into the noise timestream and almost none in the data. The bad: weighting is based on the noise timestream so it's possible that the weights aren't quite right as a result The weird: the data PSD. What's up with that? Apparently I'm preferentially subtracting certain scales but I don't know why, unless deconvolution is at fault. Edit/Update: The deconvolution is definitely at fault. Here's the same scan done without deconvolution:
It should have been obvious; the cyan in the first plot is the PSD of the deconvolution straight up, and that should have no high-frequency structure...
Paragraph for the methods paper
The raw data from Bolocam contains noise components from the atmosphere and instrument in addition to the astrophysical signal. To remove the atmospheric noise, an iterative approach was required.
- The median across all bolometers is subtracted
- A set number of principle components are subtracted. The principle components are the most correlated components between bolometers. In this process both the atmosphere above the telescope - which is assumed to be constant across the field of view - and any large-scale astrophysical structure are removed.
- The timestream data is mapped into the plane of the sky. Data points are mapped to the nearest pixel. 7.2" pixels are used so that sampling is better than Nyquist.
- The map is deconvolved using a maximum entropy deconvolution algorithm ( Based on paper by Hollis, Dorband, Yusef-Zadeh, Ap.J. Feb.1992, written by Frank Varosi at NASA/GSFC 1992)
- The deconvolved map is returned to a timestream and subtracted from the original to yield a noise-only timestream.
- Power spectral densities are calculated for each scan in the noise timestream, and weights are calculated from these. [At the moment, the weights are actually inverse-variance]
- The deconvolved map timestream is subtracted from the raw timestream, and then steps 1-6 are repeated on that timestream to recover flux that was oversubtracted in the first iteration.
Convergence takes ???? iterations.... ??? PCA components are subtracted [default 13]...
Rewriting pipeline again
Not entirely, but I realized my code had become... clunky, at best. I'm rewriting the wrappers using structs. I HATE IDL structs, but they're necessary.... if only they would make them dynamically modifiable. Weighting worked to some degree, the mapping is pretty much done.