基于艾伦方差的频率稳定性分析

某个授时系统通过串口或网口采集时间间隔计数器、频率计数器、相位噪声分析仪设备的重要信息，用于评估和分析频率源的频率稳定度，确保测量的准确性和可靠性。

数据处理：

1. 读取保存在文件中的时间间隔计数器测量的时差数据，计算时间稳定度（用TDEV表示）并保存。TDEV包括秒稳，十秒稳，百秒稳，千秒稳和万秒稳的计算。
2. 时间稳定度的计算方法：Allan方差法：将时间间隔计数器采集的数据根据艾伦方差公式计算时间稳定度。

Allan 方差反应了相邻两个采样段内平均频率差的起伏，它的最大优点在于对各类噪声的幂律谱项都是收敛的（对于那些幂律谱噪声，Allan 方差的计算结果不会出现无穷大或者无规律的发散情况）。

Allan方差是时频分析和惯性导航领域常用的一种误差分析方法，它有效地刻画了待研究的误差时间序列在不同时间尺度上的波动水平（不稳定性），并可根据不同时间尺度上的Allan方差值所构成的曲线的形状特征来辨识其中包含的随机过程模型。Allan方差分析方法对中长期的随机波动具有很强的表现力，它完全可以作为一个通用时间序列工具来推广到其它应用领域，就像PSD这样的频域分析方法和自相关这样的时域分析方法一样。

首先先理解艾伦方差的标准式分块和交叠式分块的公式，在这里贴两张自己在实现代码时找到的资料。

标准式分块计算：

交叠式分块计算：

其次是代码实现:

main.c

#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include"queue.h"
#include"allan.h"

int main()
{
    FILE *fp = fopen("1.txt","r")
    if(NULL == fp)
    {   
        perror("fail open:");
        return -1;
    }

    char sdata[256];
    char *psdata = NULL;
    char *pdata = NULL;
    double data = 0;
    
    Queue *queue = create_queue(); //创建一个队列用来储存数据
    if(NULL == queue)
    {
        printf("fail create queue");
        return -1;
    }
    
    while(1)
    {
        pdata = fgets(sdata,sizeof(sdata),fp);
        if(NULL == pdata)
        {
            break;
        }
        
        strtok(sdata," ");
        strtok(NULL," ");
        psdata = strtok(NULL," ");
        data = atof(psdata);

        push_queue(queue.data); //入队
    }

    return 0;    
}

allan.c

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include"queue.h"

int one_steady_block(Queue *queue) //计算秒稳的块数
{
    return all_datas_count(queue); //返回的函数值为数据的总数，count
}
int ten_steady_block(Queue *queue) //10s稳
{
    return all_datas_count(queue) / 10; 
}
int hundren_steady_block(Queue *queue)
{
    return all_datas_count(queue) / 100;
}
int thousand_steady_block(Queue *queue)
{
    return all_datas_count(queue) / 1000;
}
int ten_thousand_steady_block(Queue *queue)
{
    return all_datas_count(queue) / 10000;
}


double one_steady(Queue *queue) // 秒稳
{
    double y2_y1 = 0;                            //定义一个变量记录 y2 - y1 平方求和
    Que_Node *p = queue -> phead -> pnext;       //p指向第二个数据
    double y1 = queue -> phead -> data;          //y1 = 第一个数据
    double y2 = p -> data;                       //y2 = 第二个数据
    int N = one_steady_block(queue);             //N是数据块数
    int i = 0;

    for(i = 1;i <= (N - 1);++i)                  //因为第一块和第二块在循环外拿到了，所以计 
                                                   算时i从1开始，到块数 - 1
    {
        y2_y1 += (y2 - y1) * (y2 - y1);          
        p = p -> pnext;
        y1 = y2;  
        if(NULL == p)
        {
            break;
        } 
        y2 = p -> data;
    }

    return sqrt(y2_y1 / (2 * (N - 1)));
}

double ten_steady(Queue *queue) // 标准式分块的10s稳,后面的100s稳等实现方式与这个一模一样， 
                                   所以只写一个
{
    double mid = 0;                         //mid是计算y1(前位块的平均值)的一个中间值
    double y2_y1 = 0;
    double y1 = 0;
    double y2 = 0;
    int N = ten_steady_block(queue);
    Que_Node *p = queue -> phead;

    int i = 0;
    int j = 0;
    for(i = 0;i < 10;++i)                   //先在循环外拿到第一块数据的平均值
    {
        mid += p -> data;
        p = p -> pnext;
    }
    y1 = mid / 10;
    mid = 0;

    for(i = 1;i <= (N - 1);++i)
    {
        for(j = 0;j < 10;++j)               //拿到后位置上的块的平均值
        {
            mid += p -> data;
            p = p -> pnext;
        }
    }
    y2 = mid / 10;
    mid = 0;
    
    y2_y1 += (y2 -y1) * (y2 - y1);       

    y1 = y2;                                 //把后位置块的值赋给前位置块，然后再次循环
}


double ten_cross_allan(Queue *queue)
{
    double y2_y1 = 0;
    double mid = 0;
    Que_Node *p = queue -> phead;
    Que_Node *q = queue -> phead -> pnext;
    double *arr = (double *)malloc(sizeof(double) * (44200 - 10 + 1)); //一共有总数据数量 
                                                                         - 块中数据的个数 
                                                                         + 1组，所以开辟这 
                                                                         么大的空间
    int i = 0;
    int j = 0;
    for(i = 0;i < 10;++i)                        //拿到第一块的平均值放入arr【0】中
    {
        mid += p -> data;
        p = p -> pnext;
    }
    arr[0] = mid / 10.00;
    p = q;                                       //p = 用q记录的地址
    q = p -> pnext;                              //q指向p的下一个
    
    for(i = 1;i < 44200 - 10 + 1;++i)            //循环数据量 - 块中数据个数 + 1次
    {
        mid = 0;
        for(j = 0;j < 10;++j)                    //给arr【i】放入上一次p指向的下一个数据块 
                                                   的平均值 即上一次记录的q              
        {
            mid += p -> data;
            p = p -> pnext;
        }
        arr[i] = mid / 10.00;
        p = q;
        q = p -> pnext;
    }
    
    for(i = 0;i < 44200 - 10 + 1;++i)             //计算总和
    {
        y2_y1 += (arr[i + 10] - arr[i]) * (arr[i + 10] - arr[i]);
    }
         
    free(arr);

    return sqrt(y2_y1 / (2 * (44200 - 2 * 10 + 1)));
}

allan.h

#ifndef _ALLAN_H
#define _ALLAN_H

#include"queue.h"

extern .....
extern ....


#end if