谱宽参数--程序备份·

import numpy as np
import statsmodels.tsa.api as smt
import math
from sympy import *
from scipy.fft import fft
import matplotlib.pyplot as plt

def specwid(data, H):
    """
    计算谱宽参数
    :param data:数据
    :param H:上限频率
    :param fs:采样频率
    :return:e谱宽参数
    """
    N = len(data)  # 信号序列的长度
    R = smt.stattools.acf(data, nlags=N - 1, fft=True)  # 计算信号序列的自相关函数
    R = np.array(R)

    sw = abs(fft(R))
    sw = sw ** 2 * 2

    result=[]

    for i in range(0,5):

        inter=0

        c = int(H*N/20480+1)

        for m in range(0, c):
            inter += sw[m] *(m *(20480/N)) ** i
        result.append(inter)

    result = np.array(result)
    e = np.sqrt(1 - ((result[2]*result[2]) / (result[0]*result[4]))) # 计算谱宽参数
    return e


def _slide_window(rows, sw_width, sw_steps):
    '''
    函数功能：
    按指定窗口宽度和滑动步长实现单列数据截取
    --------------------------------------------------
    参数说明：
    rows：单个文件中的行数；
    sw_width：滑动窗口的窗口宽度；
    sw_steps：滑动窗口的滑动步长；
    '''
    start = 0
    s_num = (rows - sw_width) // sw_steps  # 计算滑动次数
    new_rows = sw_width + (sw_steps * s_num)  # 完整窗口包含的行数，丢弃少于窗口宽度的采样数据；

    while True:
        if (start + sw_width) > new_rows:  # 如果窗口结束索引超出最大索引，结束截取；
            return
        yield start, start + sw_width
        start += sw_steps


def read_files(path):
    '''
    读取信号文件
    :param path: 文件路径
    :param N: 信号长度
    :param len: 帧长
    :return: 信号数据
    '''
    data1 = []
    with open(path, 'r') as f1:
        str1 = f1.read()
        for item in str1.split():
            str1 = int(item[1:], 16)
            data1.append(str1 / 4096 * 3.3)
    data1 = np.array(data1)
    data1 = data1.reshape(-1,1)
    return data1

def pukuancanshudataset(data1,H):
    '''
    计算单个文件的每一帧的谱宽参数
    :param paths:路径集合
    :return: 谱宽参数集合U
    '''
    fs = 20480
    _test_list = []
    U1 = []
    for start, end in _slide_window(data1.shape[0], fs, fs):
        '''
        此处可以添加for循环或者其他方式，以实现处理多列数据
        '''
        _test_list.append(data1[start:end])
        '''
        此处可添加判断条件，以实现数组堆叠
    	'''
    for i in range(0,len(_test_list)):
        phi = specwid(_test_list[i], H)
        U1.append(phi)
    return U1

def autolabel(rects):
    '''
    定义函数来显示柱状上的数值
    :param rects:
    :return:
    '''
    for rect in rects:
        height = rect.get_height()
        plt.text(rect.get_x()+rect.get_width()/2.-0.2, 1.03*height, '%s' % float(round(height,2)))

def pukuandataset(paths,H):
    '''

    :param paths:路径集合
    :return: 谱宽参数集合U
    '''
    U1 = []
    for path in paths:
        U = read_files(path)
        phi = specwid(U, H)  # 计算近似熵
        #print(phi)
        U1.append(phi)
    return U1

def zhuzhuangtu(U1,U2,name):
    x = [1,2,3]
    bar_width = 0.2
    a = plt.bar(x=range(len(x)), height=U1, color='steelblue', label='without leakage', width=0.2)
    b = plt.bar(x=np.arange(len(x)) + bar_width, height=U2, color='indianred', label='with leakage', width=0.2)
    autolabel(a)
    autolabel(b)
    x = [1, 2,3]
    plt.xticks(np.arange(len(x)) + bar_width / 2, x)

    plt.ylabel('Apen')
    plt.xlabel('number of times')
    plt.legend()
    plt.savefig('第三次数据'+name+'Apen.jpg')
    plt.show()

def zhexiantu(U1,U2,name):
    '''

    :param U1:无泄露数据
    :param U2:有泄露数据
    :param name:折线图文件名
    :return:
    '''
    x_axis_data = [i for i in range(len(U1))]
    y_axis_data1 = U1
    y_axis_data2 = U2

    plt.clf()
    plt.plot(x_axis_data, y_axis_data1,label = 'without leakage')
    plt.plot(x_axis_data, y_axis_data2,label = 'with leakage')
    plt.ylabel('pukuancanshu')
    plt.xlabel('frames')
    plt.legend()
    plt.savefig('./第三次数据结果/' + name + 'pukuancanshu.jpg')

    plt.close('all')


def zhexiantu2(U1,U2, name):
    '''
    :param U1:无泄露数据
    :param U2:有泄露数据
    :param name:折线图文件名
    :return:
    '''
    x_axis_data = [0,1,2,3,4]
    print(x_axis_data)

    plt.clf()
    for j in range(len(U1)):
        plt.plot(x_axis_data, U1[j], color = 'b',label = 'without leakage')
        plt.plot(x_axis_data, U2[j], color = 'r',label = 'with leakage')

    plt.ylabel('pukuancanshu')
    plt.xlabel('frames')
    plt.legend()
    plt.savefig('./硬件修改后数据结果/' + name + 'pukuancanshu.jpg')

    plt.close('all')

def save_result(path1,path2,H_,beishu):
    '''
    存储运行结果
    '''
    U1 = []
    U2 = []
    for i in range(0, len(path1)):
            data1 = read_files(path1[i])
            U1.append(pukuancanshudataset(data1,H_))
            data2 = read_files(path2[i])
            U2.append(pukuancanshudataset(data2,H_))

            print(U1)
            print(U2)
            name = '积分上限频率为' + str(H_) + '/' + '放大倍数为 '+str(beishu)
            zhexiantu2(U1, U2, name)

if __name__ == "__main__":

    path10_1 = ['硬件修改后第一次听漏数据11-2/10倍不漏(第1次).txt','硬件修改后第一次听漏数据11-2/10倍不漏(第2次).txt','硬件修改后第一次听漏数据11-2/10倍不漏(第3次).txt']
    path10_2 = ['硬件修改后第一次听漏数据11-2/10倍有漏(第1次).txt','硬件修改后第一次听漏数据11-2/10倍有漏(第2次).txt','硬件修改后第一次听漏数据11-2/10倍有漏(第3次) .txt']

    path20_1 = ['硬件修改后第一次听漏数据11-2/20倍不漏(第1次).txt','硬件修改后第一次听漏数据11-2/20倍不漏(第2次).txt','硬件修改后第一次听漏数据11-2/20倍不漏(第3次).txt']
    path20_2 = ['硬件修改后第一次听漏数据11-2/20倍有漏(第1次).txt','硬件修改后第一次听漏数据11-2/20倍有漏(第2次).txt','硬件修改后第一次听漏数据11-2/20倍有漏(第3次).txt']

    path50_1 = ['硬件修改后第一次听漏数据11-2/50倍不漏(第1次).txt','硬件修改后第一次听漏数据11-2/50倍不漏(第2次).txt','硬件修改后第一次听漏数据11-2/50倍不漏(第3次).txt']
    path50_2 = ['硬件修改后第一次听漏数据11-2/50倍有漏(第1次).txt','硬件修改后第一次听漏数据11-2/50倍有漏(第2次).txt','硬件修改后第一次听漏数据11-2/50倍有漏(第3次).txt']

    path100_1 = ['硬件修改后第一次听漏数据11-2/100倍不漏(第1次).txt','硬件修改后第一次听漏数据11-2/100倍不漏(第2次).txt','硬件修改后第一次听漏数据11-2/100倍不漏(第3次).txt']
    path100_2 = ['硬件修改后第一次听漏数据11-2/100倍有漏(第1次).txt','硬件修改后第一次听漏数据11-2/100倍有漏(第2次).txt','硬件修改后第一次听漏数据11-2/100倍有漏(第3次).txt']

    H = [2000, 4000, 6000, 8000]
    for H_ in H:
        save_result(path10_1, path10_2, H_, 10)
        save_result(path20_1, path20_2, H_, 20)
        save_result(path50_1, path50_2, H_, 50)
        save_result(path100_1, path100_2, H_, 100)