# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
import pywt


#平滑噪声—等深分箱—均值平滑
def aequilatus_box_mean(data,bins):
    length=data.shape[0]
    labels=[]
    for i in range(bins):
        labels.append('a'+str(i+1))#添加标签
    new_data=pd.qcut(data.iloc[:,0],bins,labels=labels)#等深分箱
    data['label']=new_data
    for label in labels:
        label_index_min=data[data.label==label].index.min()#分箱后索引最小值
        label_index_max=data[data.label==label].index.max()#分箱后索引最大值
        data.loc[label_index_min:label_index_max,data.columns[0]]=np.mean(
            data.A[label_index_min:label_index_max+1,])#根据label及索引，修改A为各箱均值
    return data

def aequilatus_box_median(data,bins):
    length=data.shape[0]
    labels=[]
    for i in range(bins):
        labels.append('a'+str(i+1))
    new_data=pd.qcut(data.A,bins,labels=labels)#等深分箱
    data['label']=new_data
    for label in labels:
        label_index_min=data[data.label==label].index.min()#分箱后索引最小值
        label_index_max=data[data.label==label].index.max()#分箱后索引最大值
        data.loc[label_index_min:label_index_max,'A']=np.median(
            data.A[label_index_min:label_index_max+1,])#根据label及索引，修改A为各箱均值
    return data

#平滑噪声—等深分箱—边界平滑
def aequilatus_box_border(data,bins):
    length=data.shape[0]
    labels=[]
    for i in range(bins):
        labels.append('a'+str(i+1))
    new_data=pd.qcut(data.A,bins,labels=labels)#等深分箱
    data['label']=new_data
    for label in labels:
        label_index_min=data[data.label==label].index.min()
        label_index_max=data[data.label==label].index.max()
        data_min=np.min(data.A[label_index_min:label_index_max+1,])
        data_max=np.max(data.A[label_index_min:label_index_max+1,])
        for i in range(label_index_min,label_index_max):
            if(data.loc[i,'A']==data_min or data.loc[i,'A']==data_max):
                data.loc[i,'A']=data.loc[i,'A']
            elif(np.abs(data.loc[i,'A']-data_min)<=np.abs(data.loc[i,'A']-data_max)):
                 data.loc[i,'A']=data_min
            else:
                data.loc[i,'A']=data_max
    return data
#一维数据小波阈值去噪
#封装成函数
def sgn(num):
    if(num > 0.0):
        return 1.0
    elif(num == 0.0):
        return 0.0
    else:
        return -1.0

def wavelet_noising(new_df):
    data = new_df
    data = data.values.T.tolist()  # 将np.ndarray()转为列表
    w = pywt.Wavelet('sym8')#选择sym8小波基
    [ca5, cd5, cd4, cd3, cd2, cd1] = pywt.wavedec(data, w, level=5)  # 5层小波分解

    length1 = len(cd1)
    length0 = len(data)

    Cd1 = np.array(cd1)
    abs_cd1 = np.abs(Cd1)
    median_cd1 = np.median(abs_cd1)

    sigma = (1.0 / 0.6745) * median_cd1
    lamda = sigma * math.sqrt(2.0 * math.log(float(length0 ), math.e))#固定阈值计算
    usecoeffs = []
    usecoeffs.append(ca5)  # 向列表末尾添加对象

    #软硬阈值折中的方法
    a = 0.5

    for k in range(length1):
        if (abs(cd1[k]) >= lamda):
            cd1[k] = sgn(cd1[k]) * (abs(cd1[k]) - a * lamda)
        else:
            cd1[k] = 0.0

    length2 = len(cd2)
    for k in range(length2):
        if (abs(cd2[k]) >= lamda):
            cd2[k] = sgn(cd2[k]) * (abs(cd2[k]) - a * lamda)
        else:
            cd2[k] = 0.0

    length3 = len(cd3)
    for k in range(length3):
        if (abs(cd3[k]) >= lamda):
            cd3[k] = sgn(cd3[k]) * (abs(cd3[k]) - a * lamda)
        else:
            cd3[k] = 0.0

    length4 = len(cd4)
    for k in range(length4):
        if (abs(cd4[k]) >= lamda):
            cd4[k] = sgn(cd4[k]) * (abs(cd4[k]) - a * lamda)
        else:
            cd4[k] = 0.0

    length5 = len(cd5)
    for k in range(length5):
        if (abs(cd5[k]) >= lamda):
            cd5[k] = sgn(cd5[k]) * (abs(cd5[k]) - a * lamda)
        else:
            cd5[k] = 0.0

    usecoeffs.append(cd5)
    usecoeffs.append(cd4)
    usecoeffs.append(cd3)
    usecoeffs.append(cd2)
    usecoeffs.append(cd1)
    recoeffs = pywt.waverec(usecoeffs, w)#信号重构
    return recoeffs


def wavelet_noising2(new_df):
    data = new_df
    data = data.values.T.tolist()  # 将np.ndarray()转为列表
    w = pywt.Wavelet('dB10')#选择dB10小波基
    ca3, cd3, cd2, cd1 = pywt.wavedec(data, w, level=3)  # 3层小波分解
    ca3=ca3.squeeze(axis=0) #ndarray数组减维：(1，a)->(a,)
    cd3 = cd3.squeeze(axis=0)
    cd2 = cd2.squeeze(axis=0)
    cd1 = cd1.squeeze(axis=0)
    length1 = len(cd1)
    length0 = len(data[0])

    abs_cd1 = np.abs(np.array(cd1))
    median_cd1 = np.median(abs_cd1)

    sigma = (1.0 / 0.6745) * median_cd1
    lamda = sigma * math.sqrt(2.0 * math.log(float(length0 ), math.e))
    usecoeffs = []
    usecoeffs.append(ca3)

    #软阈值方法
    for k in range(length1):
        if (abs(cd1[k]) >= lamda/np.log2(2)):
            cd1[k] = sgn(cd1[k]) * (abs(cd1[k]) - lamda/np.log2(2))
        else:
            cd1[k] = 0.0

    length2 = len(cd2)
    for k in range(length2):
        if (abs(cd2[k]) >= lamda/np.log2(3)):
            cd2[k] = sgn(cd2[k]) * (abs(cd2[k]) - lamda/np.log2(3))
        else:
            cd2[k] = 0.0

    length3 = len(cd3)
    for k in range(length3):
        if (abs(cd3[k]) >= lamda/np.log2(4)):
            cd3[k] = sgn(cd3[k]) * (abs(cd3[k]) - lamda/np.log2(4))
        else:
            cd3[k] = 0.0

    usecoeffs.append(cd3)
    usecoeffs.append(cd2)
    usecoeffs.append(cd1)
    recoeffs = pywt.waverec(usecoeffs, w)#信号重构
    return recoeffs


if __name__=="__main__":
    data=pd.DataFrame({'A':[11,13,15,20,20,23,26,29,35]})
    bins=3
    print("均值平滑")
    print(aequilatus_box_mean(data,3))

    #############################################################
    print("中值平滑")
    print(aequilatus_box_median(data, 3))

    ##############################
    print("边界平滑")
    print(aequilatus_box_border(data, 3))

    #参考：   https://zhuanlan.zhihu.com/p/157540476
    #####################################################
    print("一维数据小波阈值去噪")
    data_denoising = wavelet_noising(data["A"])  # 调用函数进行小波阈值去噪
    print(data_denoising)

    #####################################################
    print("一维数据小波阈值去噪2")
    data_denoising2 = wavelet_noising2(data["A"])  # 调用函数进行小波阈值去噪
    print(data_denoising2)