APEnt_test

近似熵主要的是从衡量时间序列复杂性的角度来度量信号中产生新模式的概率大小,产生新模式的概率越大,序列的复杂性越大,相应的近似熵也越大。本篇设计了两个需要预测的时间序列,test_data1为信噪比较高的信号,test_data2为噪声较大的信号,通过计算发现test_data1的近似熵较为真实的反应了新模式的产生,而test_data2的近似熵则无法反应新模式的产生。对test_data2通过小波滤波后再进行近似熵的计算,则能较好地反应新模式的产生

import matplotlib.pyplot as plt
#from pandas import Series
#from pandas import DataFrame
#import pandas as pd
import numpy as np
##from pandas import *
#import pandas
#import xlrd
# SampEn calculation
from numpy import arange
import math
#import numpy as np
#from scipy import interpolate
#from scipy.interpolate import spline

#近似熵
def ApEn(U, m, r):

    def _maxdist(x_i, x_j):
        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

    def _phi(m):
        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
        C = [len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) / (N - m + 1.0) for x_i in x]
        return (N - m + 1.0)**(-1) * sum(np.log(C))

    N = len(U)

    return abs(_phi(m+1) - _phi(m))

#num_pi为要产生几个π的sin数据,num_ex为异常点的个数,num_gap为段缺失数据的个数,num_bk为单个缺失值的个数
def test_data_gen(num_pi,num_ex,num_gap,num_bk):
    if (num_pi>0) :
        num_point=72*num_pi
        x=np.linspace(0,np.pi*num_pi,num_point)
        signal1=[(2*math.sin(i)+1) for i in x] #产生测试用的num_pi个sin数据
        noise=np.random.normal(0,0.05,num_point)#numpy.random.normal(噪声均值, 噪声标准差, 噪声的shape)
        signal1=signal1+noise#在sin数据上添加噪声
    else:
        print("Please input valid num_pi")
        return

    if (num_ex>0) :
        #随机添加异常值
        point_ex=[]
        for i in range(num_ex):
            point_ex.append(np.random.randint(0,len(signal1))) #异常值的位置
        for _ in point_ex:
            signal1[_]=signal1[_]*1.8
    else:
        pass
    if (num_gap>0) :    
        #随机添加段数据缺失
        longth_gap=np.random.randint(15)+5 #缺口大小5~20

        point_gap=[]   #缺口的位置
        for i in range(num_gap):
            point_gap.append(np.random.randint(num_point-20))

        for i in point_gap:
            for j in range(longth_gap):  
                signal1[i+j]=None
    else:
        pass
    if (num_bk>0) :        
        #随机添加单点缺失值
        point_break=[]
        for i in range(num_bk):
            point_break.append(np.random.randint(num_point))        
        for _ in point_break:
            signal1[_]=None
    else:
        pass
    #产生时间序列,每隔5分钟一个点
    """
    date_need=[]
    start_dt = datetime.datetime(2017, 1, 1) 
    interval = datetime.timedelta(seconds=300) 
    for i in range(num_point): 
        date_need.append(start_dt + interval * i)

    df = DataFrame(signal1,index = date_need[0:num_point])
    df.to_excel('data_test.xlsx')
    """
    plt.figure(figsize=(10,5))

    plt.plot(signal1)
    plt.show()

    return signal1
test_data=test_data_gen(240,0,0,0)

png这里写图片描述

num_point=72*240
x=np.linspace(0,np.pi*3,num_point)
signal_modul=[(math.sin(i)+0.4) for i in x] #调制用的信号
for i in range(72*240):
    if (i<1000 or i>3000):
        signal_modul[i]=1
signal_modul[9000:12000]=signal_modul[1000:3000]
signal_modul[1000:3001]=[1]*3000
signal_modul.append(1)
plt.plot(signal_modul)
[<matplotlib.lines.Line2D at 0x13c753aad68>]

png这里写图片描述

data=signal_modul*test_data
plt.plot(data)
[<matplotlib.lines.Line2D at 0x13c753f2358>]

png这里写图片描述

z_list=[]
x_list=[]
j = (len(data)-288*7)//288+1
for i in range(j):

    y = data[i*288:i*288+288*7]
    y_std = np.std(y, ddof=1)#ddof=1为样本标准偏差,而非整体标准偏差
    z = ApEn(y, 2, 0.1*y_std)#将数据切片为7天一组,计算一组的近似熵
#    z = entropy(y)

    #print(i)
    #print(z)
#    print(y)
    z_list.append(z)
    x_list.append(i*288+288*7)
    x1 = arange(0, len(data), 1)
#x2 = cols0[1400:1400+j]
#fig=plt.figure(figsize=(20,10))
#fig, axes = plt.subplots(2, 1,  sharex=True,figsize=(25,10))

#ax1=fig.add_subplot(2,1,1)
#ax2=fig.add_subplot(2,1,2)
fig, axes = plt.subplots(2, 1,  sharex=True,figsize=(25,10))
axes[0].scatter(x_list,z_list)
axes[1].plot(x1,data)
[<matplotlib.lines.Line2D at 0x13c76483668>]

png这里写图片描述

一开始没有设置好test_data中噪声的幅度,导致近似熵无法很好的反应新模式的产生

#num_pi为要产生几个π的sin数据,num_ex为异常点的个数,num_gap为段缺失数据的个数,num_bk为单个缺失值的个数
def test_data_gen2(num_pi,num_ex):
    if (num_pi>0) :
        num_point=72*num_pi
        x=np.linspace(0,np.pi*num_pi,num_point)
        signal1=[(math.sin(i)+1) for i in x] #产生测试用的num_pi个sin数据
        noise=np.random.normal(0,0.1,num_point)#numpy.random.normal(噪声均值, 噪声标准差, 噪声的shape)
        signal1=signal1+noise#在sin数据上添加噪声
    else:
        print("Please input valid num_pi")
        return

    if (num_ex>0) :
        #随机添加异常值
        point_ex=[]
        for i in range(num_ex):
            point_ex.append(np.random.randint(0,len(signal1))) #异常值的位置
        for _ in point_ex:
            signal1[_]=signal1[_]*1.8
    else:
        pass
    plt.figure(figsize=(10,5))

    plt.plot(signal1)
    plt.show()
    return signal1
test_data2=test_data_gen2(240,0)

png这里写图片描述

data2=signal_modul*test_data2
z2_list=[]
x2_list=[]
j = (len(data2)-288*7)//288+1
for i in range(j):

    y = data2[i*288:i*288+288*7]
    y_std = np.std(y, ddof=1)#ddof=1为样本标准偏差,而非整体标准偏差
    z = ApEn(y, 2, 0.1*y_std)#将数据切片为7天一组,计算一组的近似熵
#    z = entropy(y)

    #print(i)
    #print(z)
#    print(y)
    z2_list.append(z)
    x2_list.append(i*288+288*7)
    x2 = arange(0, len(data2), 1)
fig, axes = plt.subplots(2, 1,  sharex=True,figsize=(25,10))
axes[0].scatter(x2_list,z2_list)
axes[1].plot(x2,data2)
[<matplotlib.lines.Line2D at 0x13c76a186d8>]

png这里写图片描述

import pywt

# 小波滤噪
def wavelet_denoising(data):
    # 小波函数取db4
    db4 = pywt.Wavelet('db4')
    #if type(data) is not types.NoneType:
    # 分解
    coeffs = pywt.wavedec(data, db4)
    # 高频系数置零
    coeffs[len(coeffs)-1] *= 0
    coeffs[len(coeffs)-2] *= 0
    coeffs[-3]*=0
    coeffs[-4]*=0

    # 重构
    meta = pywt.waverec(coeffs, db4)
    return meta
data3=wavelet_denoising(data2)
z3_list=[]
x3_list=[]

k = (len(data3)-288*7)//288+1
for i in range(k):

    y = data3[i*288:i*288+288*7]
    y_std = np.std(y, ddof=1)#ddof=1为样本标准偏差,而非整体标准偏差
    z = ApEn(y, 2, 0.1*y_std)#将数据切片为7天一组,计算一组的近似熵
#    z = entropy(y)

    #print(i)
    #print(z)
#    print(y)
    z3_list.append(z)
    x3_list.append(i*288+288*7)
    x3 = arange(0, len(data3), 1)
fig, axes = plt.subplots(2, 1,  sharex=True,figsize=(25,10))
axes[0].scatter(x3_list,z3_list)
axes[1].plot(x3,data3)
[<matplotlib.lines.Line2D at 0x13c005d5828>]

png这里写图片描述

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值