dec05 pulsar 9 level folding problem with the wapps

12aug05

Problem:

Project P1979 (snezana and joel) were doing HI absorption on pulsars using 9 level sampling and folding mode with the wapps dumping acfs. They reported a problem where a pulsar that was known to have no absorption showed large absorption with 9 level sampling and no absorption with 3 level sampling. It looked like the level correction was the culprit.

Background:

The steps in going from acf to spectrum in pulsar mode using the wapps are:

accumulate the acf in the correlator chip then read it out. For pulsars this is typically every 64 to 1024 usecs.
Combine multiple chips if 9 level or interleaved sampling is used.

this is done in one of the wapp demux routines.

correct for any bias/offsets of the 3/9 level correlator chips (e.g. -1,0,1 is recorded as 0,1,2 in the 3 level chips).
If pulsar folding mode then accumulate the acf in the correct phase bin.
Compute the digitizerThreshold/rmsVoltage using the zero lag.
3 or 9 level correct the acf
fourier transform the acf to get the spectrum (after symmeterizing or zero extending the acf).

Some notes:

If you are storing acf's, then the data is written to disc after step 4. If spectra are stored then they are written to disc after step 7.
For 9 level sampling, four 3 level correlator chips are combined to measure 9 levels. The input signal is split into two ternary digits HL (base 3). The 4 chips do the La*Lb, Ha*Hb, Ha*Lb, Hb*Lb where a is the immediate data and b is the delayed data that goes into each lag multiply. The multiplies are scaled by : 1*1, 3*3, 3*1, 1*3 corresponding to the ternary digit they belong to. This ends up scaling the acf by : 1+9+3+3=16.

level correction history:

The interim correlator 3 level correction comes from kulkarni and heiles (astronomical journal)
The interim correlator 9 level correction came from murray lewis. Only the 0 lag (total power) correction is done. The other lags are not modified.
The wapp original 3 level correction came from andy dowd and the 9 level correction was that of murray lewis.
Around 2002 fred schwab (nrao) published an algorithm for multi level correction. In the paper he included routines for 3 and 9 level correction. The wapps eventually switched to this code.
Idl routines (written by me) are able to read pulsar wapp acf data and do the kulkarni/heiles 3 level correction and murray lewis's 9 level correction (actually the 9 level idl code was written to find this problem and is probably not in the production library yet..).

Plots of the problem (and solution).

The pulsar 0950+08 was observed for about 600 seconds. Acf's were folded and written out every 20 seconds (giving 30 sets of 100 phase bins). Wapp1 was configured for 9 level sampling and wapp 3 was configured for 3 level sampling. I analyzed this data in 7 different modes. For each "mode", I averaged the 30 folded sets of phase bins. I then found the range of pulsar on, off bins and computed (pulsOn(frq)-pulsoff(frq))/median(pulsOn(frq)-pulsoff(frq)) to see the absorption. The divide scales the values to the pulsar continuum emission. Note that there is no divide by the bandpass so the spectra for the 9 level and 3 level results will be a little different (since they have different IF 50 Mhz filters). The 7 modes of analyzing the data were:

No 9 level spectral correction (i did do the 9 level total power correction)
Use foldfft with the 9 level bug in it
Use idl and murray lewis's 9 level correction to process the data.
Use foldfft with the fixed code for 9 level.
Use 3 level spectra with no correction (i did do the 3 level total power correction).
Use 3 level spectra with foldfft
Use idl and kulkarni/heiles 3 level correction.

The plots show the results of processing the data with the 7 different methods (.ps) (.pdf) :

Fig 1 Top: This shows the absorption spectrum for the 7 methods. The horizontal scale is in units of pulsar continuum. An offset has been added for plotting. The colored lines are:

black : no 9 level spectral correction.
red : foldfft with the 9 level bug.
green: foldfft with the 9 level fix.
darkBlue: 9 level correction of idl and murray lewis.
light blue: 3 level no spectral correction.
purple: 3 level correction of kulkarni/heiles.

Fig 1 bottom: Spectra from the first phase bin. The blue (9lvl_idl_bml) and black (9lvl_nocorrection) lie on top of each other. The 3 level and 9 level spectra can have different power levels since they went through different filters/amplifiers. The edges of each spectra should hit .3 (relative to the median value). This is the expected alias from the digital filters. The red plot (foldfft with bug) has spectral edges that fall below 0.
Fig 2 top: This shows the difference between 3 level spectra: foldfft (fred schwab) - (kulkarni/heiles). The difference is around 1e-7 which is close to the resolution of the floating point numbers. So these two routines are giving the same results.
Fig 2 middle: This shows 9 level spectra: foldfftCorrected - murrayLewis. There is an offset of about 3 % Tsys and a spectral variation of .0002 Tsys. If things were perfect, this difference should be the same range as the 3 level differences.
Fig 2 bottom: This shows the 9 level absorption spectrum difference (it's units are pulsar continuum rather than correlator counts). The bump that was seen in the difference is not evident. Since we are taking (onPulse-offPulse)/median(onpulse-offpulse) we are removing any offsets and linear scaling terms. So to within and additive and multiplicative constant the values are pretty close (although it would be nice to know why they aren't exactly the same).

So the two methods of doing the 3 level correction are identical. The 9 level correction are the same to within an additive and multiplicative constant. It turns out that these methods use slightly different values for the optimum digitizerThreshold/rmsvoltage. For this comparison I used the values from fred schwab. Maybe if murray redid his fit with the new optimum alfa/sigma the correction would be a little closer. You can see that the 9 level spectral correction does not make a big difference (since 9 levels is pretty close to multi bit). Only at high correlation values (those close to the 0 lag) does the correction do much.

What is causing the problem:

Snez was running in 9 level folding mode and writing uncorrected acf's to disc. She then ran the program foldfft (/share/wappsrc/utility/foldfft.c). This routine read in the folded data, called fred schwab's routines to 9 level correct the data and then fft it.

One of these routines in the level correction: thresh() (in /share/wappsrc/lib/vanvleck.c) is used to compute the digitizer threshold/rmsVoltage. It expects the 0 lag to be input without the factor of 16 removed. The online folding routine had removed this factor (see /share/wappsrc/wapp/fold.c line 109:

if( c->mode->level == NINELEV )

c->fold_count /= 16.0;

The thresh() thought that there was 16 times less power than actually present. This caused an incorrect threshold to be used in the later 9 level correction which made a big difference in the corrected spectrum.

The routine thresh() is called from: vanvleck_z(), pow9lev(), attndb9lev() all in /share/wappsrc/lib/vanvleck.c . These routines are then called from a bunch of other places.

How to fix the problem(s):

For the offline program /share/wappsrc/utility/foldfft():

I've created a version of foldfft() that works for 9 level sampling. All you need to do is multiply the zero lag passed into vanvleck_z() and pow9lev() by 16 (but only do this for 9 level sampling).

If you did 9 level online folding and you stored spectra:

Your spectra have this problem. To fix it you would have to go back to the acf domain, undo the "wrong" 9 level correction, and then redo the correct one (i'm not sure how easy this is..).

To fix online folding and recording of spectra code you need to :

div=1./(double)o->count[b] change to:

scl=(((c->head->lagformat & ~FOLD32) == FLOATSPEC ) &&
(c->mode->level == NINELEV)) ? 16. : 1.;
div= scl/(double)o->count[b];

You should not put factors of 16 in the vanvleck.c file since these routines are also used by the non-folding modes were the correct lag0 is being input.

What is not broken:

It was the folding accumulation that divided by 16. Pulsar 9 level search and spectral line 9 level do not have this division so they should be ok (or at least not have this problem).

processing: x101/051212/snez.pro, foldfft.c, avgfile.pro, inpcmp.pro

 home_~phil