SIGPROC — 脉冲星信号处理程序

Persus & Xie

原文参照 Dunc Lorimer --Jodrell Bank

SIGPROC包是一个程序的集合，主要目标是转换和处理快速采样脉冲星数据到一个紧凑和易于使用的格式，这种格式适合离线分析信号，包括信号的搜索、时间和偏振测量应用。其主要处理的数据为filterbank格式。本文档描述如何安装和运行各个命令。并利用真实数据集和模拟数据集作为应用实例。

关于sigproc

SIGPROC是一个用于对多种快速采样脉冲星数据的初始分析进行标准化的软件包。目前公认的设备有the Wide Band Arecibo Pulsar Processor (WAPP), the Penn State Pulsar Machine (PSPM), Arecibo Observatory Fourier Transform Machine (AOFTM), Berkeley Pulsar Processors (BPP), Parkes/Jodrell 1-bit filterbanks (SCAMP) 和 filterbank at the Ooty radio telescope (OOTY).这个包应该可以帮助用户快速查看他们的数据，而不需要编写另一个例程来读取数据，也不需要担心大/小端兼容性(字节交换是自动处理的)。目前的子程序是:

filterbank – 将脉冲星原始数据转换为filterbank格式：与多极化和/或频率通道相关的n-bit数字流。

splice – 将相同时间采样的filterbank数据文件并合

fake – 产生模拟的filterbank格式数据，包含存在于高斯噪声中的周期信号，用于接下来程序的测试和校准。

decimate --将输入的filterbank数据的频率和/或时间采样加在一起，以减少时间和/或频率分辨率，如果想要快速查看脉冲星脉冲轮廓或者其他信息很有帮助。

dedisperse --将输入的filterbank数据进行消色散，将输出的时间序列写成一个或者多个去色散的sub-bands.

fold – 将输入的filterbank数据或时间序列数据以脉冲周期的形式折叠。脉冲输出为ASCII EPN (§C)或psrfits格式的数据。也可以使用相应的脚本生成多项式系数。

profile – 以ASCII格式或伪灰度图显示概要文件到标准输出

pgplotter – 将来自fold和其他SIGPROC输出的概要文件显示到PGPLOT窗口

bandpass – 将平均 bandpass写入ASCII文件

header – 读取原始数据文件，filterbank或时间序列数据，用于显示表头信息，以纯ASCII形式

reader – 读取filterbank或时间序列数据，并以直接可读的形式显示

quicklook – 用csh脚本对已知脉冲星上的全功率filterbank数据进行快速分析

monitor – 使用Tk弹出小部件监视在给定目录中运行的程序脚本

管道命令详解

学习管道之前我们先了解一下Linux的命令执行顺序

命令执行顺序控制

通常情况下，我们在终端中只能执行一条命令，然后按下回车执行，那么如何执行多条命令的？

顺序执行多条命令：command1;command2;command3;

简单的顺序指令可以通过 ; 来实现

有条件的执行多条命令： command1 && command2 || command3

&& :如果前一条命令执行成功则执行下一条命令，如果command1执行成功则执行第二条命令

|| :与&&相反，是前一条命令执行不成功时执行后边的命令

$? :储存上一条命令的返回结果

实例1.1 ，在目录下有三个文件001.txt，002.txt，003.txt

la && ls
ls && la
ls || la
la || ls

运行结果

>>> command not found: la
>>> 001.txt 002.txt 003.txtcommand not found: la
>>> 001.txt 002.txt 003.txt
>>> command not found: la001.txt 002.txt 003.txt

管道命令

管道命令是一种通信机制，通常用于进程间的通讯（也可以通过socket进行通信），它表现出来的形式是将前面每一个进程的输出（stdout）直接作为下一个进程的输入（stdin）。

管道命令 | 作为界定符号，管道命令与上面说的连续执行命令不一样。

管道命令仅能处理standard output, 对于stand error output会予以忽略。

less,more,head,tail...都可以接受standard input命令，所以他们是管道命令

ls,cp,mv并不会接受standard input的命令，所以他们就不是管道命令了

管道命令必须要能够接受来自前一个命令的数据称为standard input继续处理才行

第一个管道命令

ls -al /etc | less

通过管道将ls -al的输出作为下一命令less的输入，方便浏览:

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-NwWalqbw-1597224956423)(/Users/persusx/Documents/OneDrive/文档/管道命令流程图.jpg)]

管道命令的处理图

关于sigproc 的pipe

SIGPROC中的所有程序都是用C语言编写的，可以在Linux/UNIX 系统中以命令行的形式执行命令。使用是由标准的输入和输出流组成的，这样就可以用管道| 将各种任务的程序命令“粘合”在一起。例如，下面的命令:

filterbank B0823+26.pspm | dedisperse -d 19 -s 4 | fold -p polyco.dat > B0823+26.prf

上述命令目的是：将读取原始PSPM数据并将其分解成四个子带，然后根据存储在文件polyco.dat中的TEMPO生成的一组多项式系数将脉冲周期模数折叠起来。每个波段的折叠配置文件都以ASCII格式写入B0823+26.prf文件。

这些程序和脚本的详细描述将在本文档的剩余部分给出，其结构如下:在§2中，我们描述了如何安装SIGPROC;§3描述了所有程序使用的filterbank数据和头格式;产生真实和虚假的滤波器组数据，分别在§4和§5中描述;在§6和§7中讨论了查看报头和原始数据的程序;§8描述了数据简化任务(抽取和消散);在9中描述了折叠滤波组数据以产生脉冲轮廓;§10中分别给出了一个用于快速数据分析的脚本;版本历史和未来工作的计划(§11)。补充附录处理监控程序(§A)，使用TEMPO(§B)和EPN数据格式(§C)生成polyco.dat文件。

SIGPROC 软件安装

到目前为止，SIGPROC已经成功地安装在Solaris、Linux、HP-UX和mac上。在编写程序期间，ANSII 和C语言始终被坚持使用，以便在其他操作系统上安装。安装过程如下:

软件下载地址：http://www.jb.man.ac.uk/~drl/sigproc/sigproc.tar.gz

也可以在terminal中用wget获取，命令如下
```
 wget http://www.jb.man.ac.uk/~drl/sigproc/sigproc.tar.gz
```
解压软件包到指定目录，在terminal中运行，其中version代表sigproc版本号，需要根据实际情况修改命令：
```
gunzip -c sigproc-version.tar.gz | tar xvf -
```
tar文件的内容将分发到目录sigproc-version/中。进入这个目录并运行以下命令，输入:
```
cd sigproc-version
./bootstrap
./configure --prefix=$PSRSOFT_DIR --with-cfitsio-dir=$PSRSOFT_DIR CC=gcc CXX=g++ F77=gfortran CFLAGS=-fPIC FFLAGS=-fPIC  --with-fftw-dir=$PSRSOFT_DIR  CPPFLAGS=-I$PSRSOFT_DIR/include LDFLAGS="-L$PSRSOFT_DIR/lib -L$PGPLOT_DIR -L/usr/lib64”
make && make install
make clean
```
其中在配置环境，即./configure时，需要提前对$PSRSOFT_DIR在bash文件中进行设置，另，对于macOS 系统，需要在运行make && make install时对/sigproc-version/src/swap_bytes.c文件作出修改，修改函数int little_endian(),修改内容如下：

int little_endian() /includefile/

{

char *ostype;

if((ostype = (char *)getenv(“OSTYPE”)) == NULL )

error_message(“environment variable OSTYPE not set!”);

if (strings_equal(ostype,“linux”)) return 1;

if (strings_equal(ostype,“hpux”)) return 0;

if (strings_equal(ostype,“solaris”)) return 0;

if (strings_equal(ostype,“darwin”)) return 0;

return 0;

}

注，对于macOS 系统，在配置环境时即./configure时，可能需要对gcc, g++作出调整，需要查看电脑中所对应的gcc和g++的对应命令，并作出调整。当出现提示时，提供您希望将SIGPROC可执行文件放置在其中的目录的名称。如果在多个系统上编译，请登录到另一个系统并在这台计算机上运行相同的脚本。注意，如果在多个平台下编译，只需要一个源代码副本。

filerbank 数据格式

filterbank文件包含了表头和数据，描述其数据格式的目的是为了使得数据可以在其他程序中被使用，或者是有利于写出读取、处理filterbank数据格式的程序。

其数据格式可以描述为

HEADER_START stream_of_header_parameters HEADER_END stream_of_data_values

即：

表头开始项表头参数及对应数据描述表头结束项数据流

HEADER_START和HEADER_END字符串表示描述数据的头参数流的开始和结束。默认情况下是在数据文件的开头包含这些内容。我们认识到，有些用户不喜欢以这种方式处理header。对于这些用户，filterbank有一个-headerfile命令行选项，可以将消息头传输到一个单独的ASCII文件中。

为了便于使用，表头变量设定了关键参数。目前,这些包括:

telescope_id (int): 0=fake data; 1=Arecibo; 2=Ooty… others to be added

**machine_id **(int): 0=FAKE; 1=PSPM; 2=WAPP; 3=OOTY… others to be added

data_type (int): 1=filterbank; 2=time series… others to be added

rawdatafile (char []): the name of the original data file

source_name (char []): the name of the source being observed by the telescope

barycentric (int): equals 1 if data are barycentric or 0 otherwise

pulsarcentric (int): equals 1 if data are pulsarcentric or 0 otherwise

az_start (double): telescope azimuth at start of scan (degrees)

za_start (double): telescope zenith angle at start of scan (degrees)

src_raj (double): right ascension (J2000) of source (hhmmss.s)

src_dej (double): declination (J2000) of source (ddmmss.s)

tstart (double): time stamp (MJD) of first sample

tsamp (double): time interval between samples (s)

nbits (int): number of bits per time sample

nsamples (int): number of time samples in the data file (rarely used any more)

fch1 (double): centre frequency (MHz) of first filterbank channel

foff (double): filterbank channel bandwidth (MHz)

**FREQUENCY_START **(character): start of frequency table (see below for explanation)

fchannel (double): frequency channel value (MHz)

FREQUENCY_END (character): end of frequency table (see below for explanation)

nchans (int): number of filterbank channels

nifs (int): number of seperate IF channels

refdm (double): reference dispersion measure ( cm−3cm^{-3}cm−3 pc)

period (double): folding period (s)

转换用中文介绍表头信息可能会在理解上导致偏差，所以以下内容作为参考：

telescope_id (整数):望远镜代表的id,其中0是fake生成的模拟数据, 1=Arecibo; 2=Ooty… 其他可添加

**machine_id **(整数):接收机id， 0=FAKE; 1=PSPM; 2=WAPP; 3=OOTY… 其他可添加

data_type (整数):数据类型， 1=filterbank; 2=time series… 其他可添加

rawdatafile (字符 []): 原始数据的名字

source_name (字符 []): 望远镜所观测的源的名字

barycentric (整数): 如果数据是关于有质心的则等于1，否则等于0

pulsarcentric (整数): 如果波束是以脉冲星为中心等于1，否则等于0

az_start (浮点型数据): 望远镜开始观测时的方位角 (degrees)

za_start (浮点型数据): 望远镜开始观测时的俯仰角(degrees)

src_raj (浮点型数据): 观测源的赤经 (J2000) (hhmmss.s)

src_dej (浮点型数据): 观测源的赤纬 (J2000) (ddmmss.s)

tstart (浮点型数据): 第一次采样开始的时间 (MJD)

tsamp (浮点型数据): 采样时间间隔 (s)

nbits (整数): 每个采样的bits数

nsamples (整数): 数据文件中采样的时间数 (很少被用到)

fch1 (浮点型数据): 第一个滤波通道的中心频率(MHz)

foff (浮点型数据): 滤波通道带宽(MHz)

**FREQUENCY_START **(字符): 频率表的起始 (see below for explanation)

fchannel (浮点型数据): 频率通道数值 (MHz)

FREQUENCY_END (字符): 频率表的终止(see below for explanation)

nchans (整数): 滤波通道的数目

nifs (整数): 分离的IF通道的数量

refdm (浮点型数据): 参考色散量 ( cm−3cm^{-3}cm−3 pc)

period (浮点型数据): 折叠周期 (s)

一个给定的表头将包含上面的大多数变量，但不一定会包含全部的变量信息，比如在很多的观测数据中并不会存在望远镜的初始方位角和俯仰角。

一般情况下，数据由nchans频率通道的nifs极化通道和nbit数nchans频率通道组成。随后的数据流可以被认为是N个元素的一维数组，索引在0到N−1之间运行，其中

N=nifs×nchans×nsamplesN=nifs\times nchans \times nsamples N=nifs×nchans×nsamples

nsamples是观察时长除以tsamp。因此，对于给定的IF通道i = (0,1,2,3)和频率通道c =(0…nchans−1)，样本**s =(0,1,2 . .)**的数组索引为

s×nifs×nchans+i×nchans+cs\times nifs \times nchans + i \times nchans +c s×nifs×nchans+i×nchans+c

频率通道为c的观测频率可用下式简单计算：

fch1+c×fofffch1 + c \times fofffch1+c×foff

filterbank数据遵循Parkes/Jodrell Bank约定，在表头中为foff分配一个负频率，以表示最高频率通道为fch1。目前，所有的filterbank数据都是按照这个顺序写出来的，而消色散程序在它的消色散算法中也依赖于这个情况。

虽然这个系统在大多数应用程序中使用得很好，但是从2.3版本开始，就有了一种更灵活的方法来描述频率通道。不再写入fch1和foff，现在可以按照以下方式将单个频率通道频率直接写入表头:

FREQUENCY_START f1 f2 f3 f4 FREQUENCY_END

在f1,f2……f1, f2……f1,f2……为频率通道值，单位为MHz。只要f1f1f1是最高频率(同样，这是根据消色散算法规定的)，这些频率可以是任意顺序的。splice程序使用这种频率表方法来处理下面描述的非连续数据。

ASCII 表头

正如在上文提到的，filterbank在写入数据之前，会建立一个表头。以方便SIGPROC的子程序使用这个表头中所包含的信息来处理数据。要在分析其他程序时使用它，可以调用函数read_header并链接文件read_header.c中包含的其他例程。对于那些不喜欢被这些例程困扰的人，可以在调用filterbank时使用-headerfile选项。例如:

filterbank B0823+26.pspm -headerfile > B0823+26.fil

将创建文件B0823+26.fil，在ASCII文件头中只包含filterbank的数据相关的参数名称以及相关参数。比如:

Original PSPM file: B0823+26.pspm
Sample time (us): 80.000002
Time stamp (MJD): 51740.882986111108
Number of samples/record: 512
Center freq (MHz): 430.000000
Channel band (kHz): 62.000000
Number of channels/record: 128

然后让用户在他觉得合适的时候解析这个文件。另一种获取表头信息的方法是使用下面的例子中的header程序:

filterbank B0823+26.pspm | header -tstart

它将向标准输出返回51740.882986111108。上文中列出的任何表头变量名都可以作为头程序的命令行选项。更多的细节在下文中给出。

用Python 读取filterbank格式数据

用Python 读取表头

读取表头信息需要用到的Python模块：struct

import struct
header_params = {"HEADER_START": 'flag',"telescope_id": 'i',"machine_id": 'i',"data_type": 'i', "rawdatafile": 'str',"source_name": 'str', "barycentric": 'i', "pulsarcentric": 'i', "az_start": 'd',  "za_start": 'd',  "src_raj": 'd',  "src_dej": 'd',  "tstart": 'd',  "tsamp": 'd',  "nbits": 'i', "signed": 'b', "nsamples": 'i', "nbeams": "i","ibeam": "i","fch1": 'd',  "foff": 'd',"FREQUENCY_START": 'flag',"fchannel": 'd',  "FREQUENCY_END": 'flag',"nchans": 'i', "nifs": 'i', "refdm": 'd',  "period": 'd',  "npuls": 'q',"nbins": 'i', "HEADER_END": 'flag'}def read_doubleval(filfile, stdout=False):dblval = struct.unpack('d', filfile.read(8))[0]if stdout:print("  double value = '%20.15f'"%dblval)return dblvaldef read_intval(filfile, stdout=False):intval = struct.unpack('i', filfile.read(4))[0]if stdout:print("  int value = '%d'"%intval)return intvaldef read_string(filfile, stdout=False):strlen = struct.unpack('i', filfile.read(4))[0]strval = filfile.read(strlen)if stdout:print("  string = '%s'"%strval)return strval.decode('utf-8')def read_paramname(filfile, stdout=False):paramname = read_string(filfile, stdout=False)print("Read '%s'"%paramname)if stdout:print("Read '%s'"%paramname)return paramnamedef read_hdr_val(filfile, stdout=False):paramname = read_paramname(filfile, stdout)print(paramname)try:if header_params[paramname] == 'd':return paramname, read_doubleval(filfile, stdout)elif header_params[paramname] == 'i':return paramname, read_intval(filfile, stdout)elif header_params[paramname] == 'q':return paramname, read_longintval(filfile, stdout)elif header_params[paramname] == 'b':return paramname, read_charval(filfile, stdout)elif header_params[paramname] == 'str':return paramname, read_string(filfile, stdout)elif header_params[paramname] == 'flag':return paramname, Noneexcept KeyError:warnings.warn("key '%s' is unknown!" % paramname)return None, Nonedef read_header(filename, verbose=False):"""Read the header of a filterbank file, and returna dictionary of header paramters and the header'ssize in bytes.Inputs:filename: Name of the filterbank file.verbose: If True, be verbose. (Default: be quiet)Outputs:header: A dictionary of header paramters.header_size: The size of the header in bytes."""header = {}filfile = open(filename, 'rb')        #for i in filfile:#   print(i)filfile.seek(0)paramname = ""while (paramname != 'HEADER_END'):if verbose:print("File location: %d" % filfile.tell())paramname, val = read_hdr_val(filfile, stdout=verbose)if verbose:print("Read param %s (value: %s)" % (paramname, val))if paramname not in ["HEADER_START", "HEADER_END"]:header[paramname] = valheader_size = filfile.tell()filfile.close()return header, header_size
header, headersize = read_header("./GBT_Lband_PSR.fil")
for i in header:print(i,header[i])

程序解析：

创建一个字典header_params，其中key为表头参数名称，value为参数的类型。

创建一个空字典header

用二进制格式读取filterbank格式文件

用seek()移动文件读取指针到指定起始位置

定义一个为字符为空的变量paramname

创建一个while循环体（当参数变量paramname不是HEADER_END时）：

先以4个字节从文件起始位置读取第一个参数名称

通过读取的参数名称和字典header_params获得参数的类型

通过获取的参数类型，再选择读取参数的值的方式并读取参数值

将获得的参数名称和参数值放入字典header

最后返回字典header

逐项打印参数名称和参数值

注：filterbank格式的数据，其表头的数据类型是：C语言结构体，所以在用Python读取的时候需要解包，即用struct模块进行解包。

用 `header`命令查看表头信息

header程序允许人们方便地访问原始数据文件，或filterbank数据格式的二进制表头字符串。作为完整输出的例子，这里是我们的GBT单脉冲测试数据的表头，另测试数据的下载地址为：GBT_Lband_PSR.fil。命令示例及结果如下：

header GBT_Lband_PSR.fil

输出结果为：

Data file                        : GBT_Lband_PSR.fil
Header size (bytes)              : 365
Data size (bytes)                : 5963411
Data type                        : filterbank (topocentric)
Telescope                        : GBT
Datataking Machine               : BPP
Source Name                      : Mystery_PSR
Source RA (J2000)                : 16:43:38.1
Source DEC (J2000)               : -12:24:58.7
Frequency of channel 1 (MHz)     : 1447.500000
Channel bandwidth      (MHz)     : -1.000000
Number of channels               : 96
Number of beams                  : 1
Beam number                      : 0
Time stamp of first sample (MJD) : 53010.484826388893
Gregorian date (YYYY/MM/DD)      : 2004/01/06
Sample time (us)                 : 72.00000
Number of samples                : 124237
Observation length (seconds)     : 8.9
Number of bits per sample        : 4
Number of IFs                    : 1

另外，header还可以与下面的一个或多个命令行选项一起使用，只返回想要的参数的值(这在从脚本中获取值而不必解析标准输出时特别有用)。目前可用的选项有:

命令选项	作用
-telescope	返回望远镜的名字
-machine	返回采样数据机器的名字
-fch1	返回第一个通道的频率单位：MHz
-foff	返回子带宽（即通道带宽）的值单位：MHz
-nchans	返回通道的数目
-tstart	返回第一个样本数据的采样时间（MJD）
-tsamp	返回采样间隔时间（us）
-nbits	返回每个样本数据的所占 bits数
-nifs	返回IF通道数目
-headersize	返回表头的大小，单位：字节（bytes）
-datasize	返回数据的大小，单位：字节（bytes）
-nsamples	返回样本数据的数量
-tobs	返回观测的时常，单位：s

需要注意的是，headersize，datasize，nsamples和tobs本身都不是header变量;它们是由程序根据文件大小和实际头变量导出的。

使用filterbank和splice进行数据转换

原始数据和SIGPROC包的其余部分之间的接口是filterbank程序。与所有程序一样，帮助信息是通过键入程序的名称，然后在后面加help:

比如：获取命令filterbank的help信息，在terminal输入命令：

filterbank help

运行结果：

filterbank - convert raw pulsar-machine data to filterbank formatusage: filterbank <rawdatafile1> .... <rawdatafileN> -{options}rawdatafile - raw data file (recognized machines: WAPP, PSPM, OOTY)options:-o filename - output file containing filterbank data (def=stdout)
-s skiptime - skip the first skiptime (s) of data (def=0.0)
-r readtime - read readtime (s) of data (def=all)
-i IFstream - write IFstream (IFstream=1,2,3,4)
-n nbits    - write n-bit numbers (def=input format)
-c minvalue - clip DM=0 samples > mean+minvalue*sigma (def=noclip)
-swapout    - perform byte swapping on output data (def=native)
-floats     - write floating-point numbers (equal to -n 32)
-sumifs     - sum IFs 1+2 to form total-power data
-headerfile - write header parameters to an ASCII file (head)
-headeronly - write ONLY binary header parametersoptions for correlator (currently WAPP) data:-hamming    - apply Hamming window before FFT (def=nowindow)
-hanning    - apply Hanning window before FFT (def=nowindow)
-novanvleck - don't do van Vleck correction before FFT (def=doit)
-invert     - invert the band after FFT (def=noinversion)
-zerolag    - write just the zero-lag value for each IF
-rawcfs     - write raw correlation functions (novanvleck)
-corcfs     - write corrected correlation functions (vanvleck)

只要给定原始数据文件的名称作为参数，filterbank就会确定数据的来源，如果它可以读取文件，在写入表头参数和数据之前解压样本，如§3所述。表头和数据为默认标准输出，但可以重定向到一个文件使用-o文件名选项，或用标准的方式:

filterbank rawdatafile > filterbankfile

由于没有其他选项，filterbank将读取和整理原始文件中的所有数据。如果只想处理指定的一部分数据，可
可以使用-r和-s命令行选项，指定数据。例如:

filterbank rawdatafile -r 10.0 > filterbankfile

只读前10秒的数据。这些选项对于快速查看数据很有帮助。

选择IF数据流或者对其进行求和

默认情况下，将读取和处理文件中的所有IF流(如果有多个)。如何选择一个或多个，忽略其他，则可以使用-i选项，示例如下：

filterbank rawdatafile -i 1 -i 2 > filterbankfile

将只处理原始数据文件的前两个IF通道。filterbank通过-sumifs选项提供了仅对前两个IF通道进行求和的选项(以形成总功率数据):

filterbank rawdatafile -sumifs > filterbankfile

这是一个有用的，例如，从偏振观测数据得到总功率，然后搜索脉冲信号。

改变每个样本的比特（bits）数

默认情况下，filterbank将以与本机格式相同的每个样本位数写入输出数据(例如PSPM的4bits采样)。对于采样bits很大的机器(例如WAPP)，使用-n选项更有效地打包数据是有用的。例如，命令:

filterbank wappdatafile -n 8 > filterbankfile

将处理WAPP数据文件(通常为每个样本16bits)，并将输出的样本打包为单字节整数。对于搜索目的，如果只看到微小的灵敏度损失，并且数据量显著减少，强烈建议使用此选项。对于WAPP数据，从16位到8位的灵敏度损失可以忽略不计，而压缩到4位则会导致约为5%的损失。

浮点型输出

目前，如果没有descaling参数是在头时，打包数据。这意味着对于需要数据的绝对值的应用程序(例如极化工作)，有必要将数据存储为浮点数。为此提供了-float选项(尽管它实际上只是-n 32的别名)。

Correlator-specific选项

目前，WAPP是SIGPROC所识别的唯一一种能够对给定的滞后数进行自动记录和偏振模式下互相关函数的相关器。自相关函数R(τ)R(\tau)R(τ)作为滞后量的函数，定义为:

R(τ)=lim⁡T→∞1T∫0+∞V(t)V∗(t+τ)dtR(\tau)=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{+\infty}V(t)V^{*}(t+\tau)dtR(τ)=T→∞limT1∫0+∞V(t)V∗(t+τ)dt

其中V (t)为复杂的信号电压对时间t的函数，由Weiner-Khinchin定理，功率谱密度P(f)P(f)P(f)为R(τ)R(\tau)R(τ)傅里叶变换后的函数:

P(f)=12π∫−∞+∞R(τ)e−2πifτdτP(f)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}R(\tau)e^{-2\pi if\tau}d\tauP(f)=2π1∫−∞+∞R(τ)e−2πifτdτ.

在实践中，获得filterbank数据的频率通道,对有限能级量子化，每个IF通道的滞后需要被纠正前的快速傅里叶变换(FFT)得到光谱——所谓的 van Vleck correction(见例如ee for example Hagen & Farley 1973, Radio Science, 8, 775–784)。作为参考，三级van Vleck公式在filterbank中用于将测量的自相关值®修正为无偏相关值(ρ\rhoρ)的可以写成:

r=1π∫0ρ(exp⁡(−(α/σ)21+x)+exp⁡(−(α/σ)2)1−x))dx(1−x2)r=\frac{1}{\pi}\int_{0}^{\rho}\left(\exp \left(\frac{-(\alpha/\sigma)^2}{1+x}\right)+\exp\left(\frac{-(\alpha/\sigma)^2)}{1-x}\right) \right)\frac{dx}{\sqrt{(1-x^2)}}r=π1∫0ρ(exp(1+x−(α/σ)2)+exp(1−x−(α/σ)2)))(1−x2)dx,

在这里α\alphaα是数字化仪的阈值，σ\sigmaσ为均方根电压。这种校正是filterbank在调用相关函数以产生光谱之前的默认做法。

有许多选项可用于修改默认处理。为了减少FFT损失，可以通过-hamming和-hanning选项命令对相关函数应用Hanning或Hamming窗口。选择-rawcfs输出量化到nbits指定的精度的原始关联函数。为了得到原始的相关函数，包括浮点数选项:

filterbank wappdatafile -rawcfs -floats > rawcffile

-corcfs选项将导出应用van Vleck校正的相关函数。

模糊相关器的选项

为了完整起见，我们提到另外两个相关器特定选项:-novanvleck和-zerolag。在FFT之前，-novanvleck选项不会应用量化校正。这个特性是真正的教学目的，因为，FFT数据得到频率通道，信噪比将丢失，如果van Vleck校正不适用。另一个主要用于测试的选项是-zerolag。如果选择，它将为每个If(所谓的零延迟)输出第一个相关函数，作为浮点数。在P(f)P(f)P(f)的上述表达式中插入τ=0\tau = 0τ=0，我们注意到零滞后只是所有频率通道的和，等价于一个没有色散测量校正的时间序列。

对于WAPP数据，最后一个选项是-invert，它将FFT后的频带反转以改变频率顺序。这通常应该在WAPP报头中处理，但在这里包含用于处理关于频率排序的报头信息不正确的数据。

拼接文件

大多数数据采集系统将每次观测收集到的数据存储为单个文件。对于阿雷西博的多个WAPP系统来说，每个系统都独立运行，对频段的不同部分进行采样，每个频段都会产生大量的数据文件。为了一起分析这些数据集，splice程序将拼合多个filterbank文件，如果他们都有相同的时间采样。命令语法非常简单:

splice file1.fil file2.fil file3.fil > splice.fil

其中假设输入文件为file1.fil, file2.fil和file3.fil，这些文件已经被转换成上面描述的filterbank格式。输出文件为splice.fil。在本例中，splice.fil也是filterbank格式，可以被后续程序读取。尽管这些文件不需要跨越一个连续的无线电频段，但是如果输入文件不具有相同的采样时间，或者它们没有按照相同频率降序排列，splice将会发报错。后一种检查是为了使数据符合消色散算法所期望的顺序。

使用fake创建模拟数据

fake程序可以用来创建包含隐藏在高斯噪声中的脉冲的测试数据集，其帮助信息为：

$ fake helpfake - produce fake filterbank format data for testing downstream codeusage: fake -{options}options:-period   p - period of fake pulsar in ms (def=random)
-width    w - pulse width in percent (def=4)
-snrpeak  s - signal-to-noise ratio of single pulse (def=1.0)
-dm       d - dispersion measure of fake pulsar (def=random)
-nbits    b - number of bits per sample (def=4)n - number of filterbank channels (def=128)
-tsamp    t - sampling time in us (def=80)
-tobs     t - observation time in s (def=10)
-tstart   t - MJD time stamp of first sample (def=50000.0)
-nifs     n - number of IFs (def=1)
-fch1     f - frequency of channel 1 in MHz (def=433.968)
-foff     f - channel bandwidth in MHz (def=0.062)
-seed     s - seed for Numerical Recipes ran1 (def=seconds since midnight)
-rednoise r - add red noise of magnitude 'r'
-nosmear    - do not add in dispersion/sampling smearing (def=add)
-swapout    - perform byte swapping on output data (def=native)
-evenodd    - even channels=1 odd channels=0 (def=noise+signal)
-headerless - do not write header info at start of file (def=header)binary options:-binary     - create binary system
-bper       - orbital period in hours (def=10.0)
-becc       - eccentricity (def=0.0, circular)
-binc       - inclination in degrees (def=90.0)
-bomega     - longitude of periastron in degrees (def= 0.0)
-bphase     - starting orbital phase (number between 0 and 1, def=0.0)
-bpmass     - pulsar mass in solar units (def=1.4))
-bcmass     - companion mass in solar units (def=5.0)

默认参数是生成一个类似于PSPM的filterbak数据。举个例子，考虑一些模拟的PSPM数据，这些数据来自于脉冲星42秒的观测，脉冲星的周期为～πms\pi\ msπ ms，占空比为10%，DM为30:

fake -period 3.1415927 -width 10 -dm 30 -tobs 42 -nbits 4 > pspm.fil

每个通道的模拟数据有一个0的平均值和单位均方根。信噪比是指每个通道中单个脉冲的高度。在上面的示例中，使用了默认的信噪比。较弱的脉冲信号可以用来实验离线搜索算法，并得出其最低限制。默认情况下，假脉冲宽度www的值依赖于filterbank设置，其中使用求和:

w2+tsamp2+tDM2\sqrt{w^2 + t_{samp}^2 +t_{DM}^2}w2+tsamp2+tDM2

其中tDMt_{DM}tDM是单个频率通道中所导致的色散展宽，可以用下式计算：

tDM=8.3×106msDMΔv/v3t_{DM}=8.3\times 10^6 ms\ DM\ \Delta v/v^3tDM=8.3×106ms DM Δv/v3,

假设通道的中心频率vvv远大于通道带宽Δv\Delta vΔv，两者单位为MHz，脉冲展宽可以用命令选项-nosnear取消掉。 Bit-format 和byte-swapping 选项和上文所描述filterbank程序中的一致。随机数生成器的起始位置默认为从午夜之后的秒数开始并多次调用随机数生成器获得的数字。这可以通过使用-seed选项指定一个位置来覆盖。

使用 `bandpass`, `reader`, `pgplotter`, `plot.py`查看数据

bandpass程序是一个简单的实用程序，读取传入数据并输出的各个通道中的均值，其均值是通过对每个通道中的所有观测时间上的值求和再平均得到的，查看bandpass命令的帮助信息：

bandpass help

输出结果如下：

bandpass - outputs the pass band from a filterbank fileusage: bandpass {filename} -{options}options:filename - filterbank data file (def=stdin)
-d numdumps - number of dumps to average over (def=all)
-t dumptime - number of seconds to average over (def=all)

在其最简单的形式，bandpass平均是对整个数据文件。图1的数据通过以下方法得到:

首先，利用fake生成周期脉冲数据：

fake -period 3.1415 -width 10 -dm 30 -tobs 120 -nifs 1 -fch1 1200 -nchans 512 -foff 0.4 -nbits 4 > pspm.fil

然后，通过bandpass程序对模拟的数据通道进行平均化处理，导出的为ASCII数据：

bandpass pspm.fil > pspm_bandpass.ascii

查看数据可以用命令pgplotter画图查看，也可以自己写python脚本查看。本案例是用python程序生成的图片。

图1 对模拟中通道数据输出

ASCII数据是写在一个简单的格式中，每一行为每个IF频率通道:if1 if2…if(nifs)。-d和-t选项允许为给定的转储数或时长计算带通的平均和输出。每个转储封装在**#START和#STOP**分隔符中:

#START
freq(1)        if(1)    if(2) ··· if(nifs)·              ·        ·          ··              ·        ·          ··              ·        ·          ·
freq(nchans)   if(1)    if(2) ··· if(nifs)
#STOP

其中**freq(1)**为第一个通道的观测频率(MHz)，以此类推，适用于所有nchans频道。SIGPROC提供了一个小的PGPLOT实用程序pgplotter，它可以绘制以这种格式传递的数据流。例如,试着

bandpass pspm.fil | pgplotter

另一个有用的程序是reader，它可以将filterbank格式的数据转换为ASCII数据流，打印到标准输出，先查看reader的帮助信息：

reader help

输出为：

reader  - look at filterbank data in ASCII formatusage: reader {filename} -{options}filename is the filterbank data file (def=stdin)options:-c        c - output only frequency channel c (1...nchans) (def=all)
-i        i - output only IF channel i (1...nifs) (def=all)
-numerate   - precede each dump with sample number (def=time)
-noindex    - do not precede each dump with number/time
-stream     - produce a stream of numbers with START/STOP boundaries

在一般情况下，一个通道数为nchans和IF通道数为nifs的filterbank文件，输出如下:

$ reader filterbankfiletime(1)  if(1)c(1)  if(1)c(2) .... if(1)c(nchans) ...... if(nifs)c(nchans)
time(2)  if(1)c(1)  if(1)c(2) .... if(1)c(nchans) ...... if(nifs)c(nchans)
time(3)  if(1)c(1)  if(1)c(2) .... if(1)c(nchans) ...... if(nifs)c(nchans)

默认情况是打印出所有IF和频率通道。输出可以通过-i和-c选项进行调整，以获得感兴趣的特定通道。例如:

reader GBT_Lband_PSR.fil -c 1 -c 2 -c 3 -c 4 | head

输出为：

0.000000 7.000000 15.000000 7.000000 8.000000
0.000072 11.000000 15.000000 9.000000 7.000000
0.000144 10.000000 15.000000 3.000000 6.000000
0.000216 9.000000 15.000000 5.000000 8.000000
0.000288 9.000000 15.000000 6.000000 4.000000
0.000360 12.000000 15.000000 6.000000 7.000000
0.000432 11.000000 15.000000 10.000000 5.000000
0.000504 12.000000 15.000000 7.000000 4.000000
0.000576 7.000000 15.000000 6.000000 11.000000
0.000648 13.000000 15.000000 8.000000 9.000000

仅显示GBT搜寻数据的前四个频率通道作为时间的函数。-numerate开关将把这个采样时间更改为一个整数计数器。时间或整数计数器可以通过-noindex选项完全关闭。-stream选项将与上面连续带通输出的情况一样，输出由#START和#STOP分隔符封装的数字流。与前面一样，可以将此格式传递给pgplotter进行绘图。

未完待续········