信噪比计算程序

如下计算信噪比的方法有两种，一种是计算峰值信噪比，一种是信号平均值信噪比。

两种都需要先对计算非脉冲信号区域的算术平方根均值，计算方法为，先定义非脉冲区域，然后将所有值求平方和，然后处以所有非脉冲区域值的个数减去一，最后求其平方根。

峰值信噪比：选取信号时间序列中的最大值处以上述所计算的非脉冲区域算术平方根均值。

信号平均值信噪比：选取信号区域，并求其均值，将求得的均值处以上述所计算的非脉冲区域算术平方根均值。

c语言

#include <string.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
int main(int argc, char *argv[])
{
  /* number of sub-integration */
  int nsub = 10;   /*  定义nsub为10   */ 
  /* number of phase bin  */
  float max[10];  /*空间为10的 max列表 */
  int mchan[10]; /* 空间为10的 mchan列表 */
  int nbin=1024 /* 定义nbin为10 */
  float peaksnr[10];  /* 空间为10的peaksnr列表, 数据类型为浮点型*/
  float sub[10][1024]; /* 空间为[10,1024]的sub列表, 数据类型为浮点型*/
  float  bin[10][1024]; /* 空间为[10,1024]的bin列表, 数据类型为浮点型*/
  float zero[10][1024]; /* 空间为[10,1024]的zero列表, 数据类型为浮点型*/
  int i,j;  /*定义整数 i,j */
  float pro[10][1024]; /* 空间为[10,1024]的pro列表, 数据类型为浮点型*/
  float   onp[10]; /* 空间为10的onp列表, 数据类型为浮点型*/
  float   offp[10]; /* 空间为10的offp列表, 数据类型为浮点型*/
  float  snr[10]; /* 空间为10的snr列表, 数据类型为浮点型*/
  /* off pulse and on pulse phase bin  */
   int onbin1,onbin2; /*定义整数 onbin1,onbin2 on_pulse区间*/
   int  offbin1=1,offbin2=800; /*定义整数 offbin1,offbin2 off_pulse区间 */
  FILE *fin; /*定义文件fin */
  FILE *fout; /*定义文件fout*/
 fin = fopen(argv[1],"r"); //读入数据文件
 printf("%s\n",argv[1]);
       
      while (!feof(fin)){ 
           //选取每个单脉冲序列中的最大值
           for (i=0; i < nsub;i++){
                max[i]=0;
              for (j=0; j < nbin;j++){ 
                 if(fscanf(fin,"%f %f %f %f",&sub[i][j],&zero[i][j],&bin[i][j],&pro[i][j])==4){
 //                  if(pro[i][j]>max[i]){
  //                   max[i]=pro[i][j];
                     //mchan[i]=j;
   //                 }
                 }
                   if(pro[i][j]>max[i]){
                     max[i]=pro[i][j];
                     mchan[i]=j;
                    }
 
             }
//               printf("%g  %d\n",max[i],mchan[i]);
           } 
       } 
   fclose(fin); 
             printf("mchan offpulse  onpulse  peak  snr   peaksnr\n");
    for (i=0; i < nsub ;i++){
           /* define the on pulse bin 定义on_pulse区间*/
           onbin1=mchan[i]-5;
           onbin2=mchan[i]+5;
       for (j=onbin1; j < onbin2;j++){ 
            onp[i]+=pro[i][j];
          } //将on_pulse区间的能量值求和
       for (j=offbin1; j < offbin2;j++){ 
                           
           offp[i]+=pro[i][j]*pro[i][j];
          } //将off_pulse区间的能量值求平方和
              onp[i]=onp[i]/(onbin2-onbin1); //对每个单脉冲时间序列on_pulse的能量求均值
              offp[i]=sqrt(offp[i]/(offbin2-offbin1-1)); //对每个单脉冲时间序列的off_pulse平方和求方根均值
              snr[i]=onp[i]/offp[i]; //	求的每个脉冲的信噪比
              peaksnr[i]=max[i]/offp[i]; //另一种峰值信噪比
              printf("%d %g %g %g %g %g\n",mchan[i],offp[i],onp[i],max[i],snr[i],peaksnr[i]);
        }
}

Python

import math
from numpy import *
import numpy as np
## 读取pulsar的单脉冲轮廓文件
data = np.genfromtxt('tutorial2.txt',dtype='float64')
## 读取 第四列：flux
a = data[:,3]
## 设定每个单脉冲的轮廓的采样点
n = 1024
## 计算脉冲的数目
m = int(len(data)/n)
print("the number of pulses is:", m)
# 设定数据维度 [脉冲数目，每个脉冲的采样点]
pulses = a.reshape(m,n)
#### 设置on pulse区间 ####
OnP_mi = 450
OnP_ma = 540
## 设置off pulse区间
off_mi1 = 0
off_ma1 = 450
off_mi2 = 541
off_ma2 = 1024
#选取每个脉冲序列中最大的值到列表pul_max
pul_max = []
for i in range(m):
    pul_max.append(max(pulses[i])

                   
pul_on = [] ##on pulse能量均值
pul_off = []  ## off pulse能量 算术平方根均值
snr = []   ## 信噪比
peaksnr = []  ##峰值信噪比
for i in range(m):
    on_e = sum(pulses[i][OnP_mi:OnP_ma])/(OnP_ma-OnP_mi)
    pul_on.append(on_e)  
    for j in np.arange(off_mi1,off_ma1):
        off_e1 += (pulses[i][j])*(pulses[i][j])
    for k in np.arange(off_mi2,off_ma2):
        off_e2 += (pulses[i][k])*(pulses[i][k])
    off_e = math.sqrt((off_e1+off_e2)/((off_ma1-off_mi1)+(off_ma2-off_mi2)-1))
    pul_off.append(off_e)
    snri = on_e/off_e
    peaksnri = pul_max[i]/off_e
    snr.append(snri)
    peaksnr.append(peaksnri)