Python【十七】：自相关函数的实现

sky

  大家可能都使用过matlab自带的xcorr函数，作为一个自相关库函数，可以对输入求其自相关序列，配合使用plot函数还可以画出图像。但是不足之处是如果想获得自相关图像中某个点的值，使用matlab是比较困难的，虽然在图像上有鼠标工具，但是如何使用程序自动实现求图像上任意一个点的值呢？一种思路是可以在相应的横坐标上画一条垂直的直线，然后求直线与图像的交点。但是这个方法对于matlab外行的自己来说，并非那么浅显明白，另外，完全封闭的xcorr函数也让自己感觉面对的是一个“黑盒子”，于是自己萌发了写一个自己使用的自相关函数的想法。由于自相关函数在不同的领域有不同的定义，我们这里实现的是一个统计学上的通用定义，即：

     这里自己使用Python来实现，需要注意，如果想针对一个文件输入求其相关性，正好适应Python将字节作为操作单位的特性，具体的代码如下：


    # -*- coding: cp936 -*-
    """
    This is a Simple Code for compute the specific value of autocorrelation function
    Compute Fomular: R(k) = E[(Xi - e)(X(i+k) - e)] / V, as E is Expection function,
    e is Expection Value and V is Variance
    e = E(Xi) = sum(Xi) / n, i = 1, 2,3 ...n
    V = E[(Xi - e)^2]
    The Line-Function of Two-point is y = Y0 + k(X - X0), k = (Y1 - Y0)/(X1 - X0)
    """

    #Autorized by GY, date 2014-08-05
    """
    本函数的作用是通过求解自相关函数在0和1处的函数值，构造其直线图像获取0附件的值，比如R(0.001)
    """
    import sys
    class File:
        def __init__(self):
            #open a file to read
            FileName = raw_input('Please input the FileName: ')
            self.fs = open(FileName, 'rb') #All attribute must be self-*
        
        def FRead(self):
            global buff
            buff = self.fs.read()
            #返回文件大小
            global size
            size = self.fs.tell()
            print 'size is ', size
            self.fs.close()
            return size

    class AR:
        #计算期望
        def Expe(self, buff):
            print 'Compute the Expection...'
            esum = 0.0
            for i in range(size):
                esum = esum + ord(buff)
            self.E = esum / size
            print 'The Expection is ', self.E
            return self.E
        #计算方差
        def Vari(self, buff):
            print 'Compute the Variance...'
            vsum = 0.0
            for e in buff:
                vsum = vsum + pow((ord(e) - self.E), 2)
            self.V = vsum / size
            print 'The Variance is ', self.V
            return self.V
        #计算坐标1点的AR值,因为根据自相关函数性质AR(0) = 1
        def AR_P(self):
            print 'compute the AR-Value...'
            self.tsum = 0.0
            for i in range(size - 1):
                self.tsum += (ord(buff)- self.E)*(ord(buff[i+1]) - self.E)
            print 'self.tsum is ', self.tsum
            self.ar = self.tsum / ((size - 1)* self.V )
            print 'self.ar is :', self.ar
            return self.ar
        #已知两点计算直线某点的y值
        def lineP(self, x):
            print 'compute the y value...'
            self.value = self.ar + (self.ar - 1)*(x - 1)
            print 'AR(',x,') is ', self.value
            return self.value

    print 'This tool is to return AutoFunction(point)...'
    fs = File()
    if fs.FRead() <= 0:
        print 'File Read Error!'
        sys.exit()
    print 'Now we begin build the AF value...'
    ar = AR()
    ar.Expe(buff)
    print ''
    ar.Vari(buff)
    ar.AR_P()
    x = input("Please input the wanted x value...")
    print 'The AR() value is '
    ar.lineP(x)
    raw_input('Enter for Exit...')


  程序其实也不长，因为毕竟只是实现了一个比较简单的自相关函数定义，下面具体来解释分析一下。
1. 第十六行到第三十行定义了一个文件类，因为这里程序的“对象”只有两个：文件输入和自相关算法，因此我们首先定义实现一个文件类，用于获得要计算相关性的文件，并且建立缓冲区存储这些输入数据；第十七行定义了一个构造函数，用于提示用户输入指定的文件名，并打开该文件获取内容；第二十二行定义了Read方法，读入内容，并且记录下文件大小（字节数），这些都会作为数据输入提供给下面的算法类，注意这里使用了全局变量global使得算法类也可以访问文件大小和输入缓冲区；
2. 第三十二行到第六十六行定义实现我们的算法类AR（AutoRelation），根据公式，我们需要分别实现期望Expe、方差Vari、AR(1)值以及直线方程四个方法。期望和方差是计算自相关函数所必须的，而直线方程式由于自相关函数是关于整数点延迟的函数，因此如果想获得0.1处的坐标只能利用直线进行估计，而这也是matlab使用的方法；
3. 第三十四行到第四十一行实现了期望函数，这里需要主要的是使用了for循环，将读入的文件的每一个字节转换成其ascii码值，利用这个值做十进制的运算加法后取平均值，为了使得结果更加精确，我们在运算中定义了浮点数；
4. 第四十二行到第五十行实现了方差函数，这里的关键是使用pow函数表示平方，然后去平均数即可；
5. 第五十二行到第六十行来计算AR(1)的值，主要思路是将x=1带入上面的公式，延迟为1即可；
6. 第六十一到第六十六行实现了直线方程，利用两点式来获取其他点的坐标，但是这不是一个好方法，对于0-1之间的点式适用的，不适用与1之后的点，那时可能需要AR(2)的值；
7. 第六十八行到第八十二行主要是程序的执行流程，并且输出一些调试信息，显示程序的运行结果；
     好了，现在来看看程序运行的结果，我们找了一个示例文件进行测试：


【PS：本节要点
1. Python按照字节为单位操作数据，非常适合字符编程；
2. 字节的ascii码值可以使用ord获得，而一个ascii码值对应的字符可以使用chr获得；
3. 面向对象的“对象”指的主体和算法；
4. 类的属性默认为局部的，即如果单独定义，则为类的所有实例共享；如果想变成全局的，必须使用global标签；
5. 引用自身的属性必须使用self前缀；
】