import pandas as pd
import numpy as np
from scipy import stats
import csv
import os
import seaborn as sns; sns.set()
#--------------------------------------时域特征提取------------------------------------------------
def get_time_domain_feature(data, window_size, step):
    """
    提取 15个 时域特征
    data, window_size, step
    @param data: shape 为 (m, n) 的 2D array 数据，其中，m 为样本个数， n 为样本（信号）长度
    @return: shape 为 (m, 15)  的 2D array 数据，其中，m 为样本个数。即 每个样本的16个时域特征
    """
    features =[]
    rows, cols = data.shape
    for i in range(0,len(data) - window_size + 1, step):
        window = data[i:i + window_size]

    # 有量纲统计量
        max_value = np.max(window)  # 最大值
        peak_value = np.max(abs(window))  # 最大绝对值
        min_value = np.min(window)  # 最小值
        mean = np.mean(window)  # 均值
    #     p_p_value = max_value - min_value  # 峰峰值
    #     abs_mean = np.mean(abs(window))  # 绝对平均值
    #     rms = np.sqrt(np.sum(window ** 2) / cols)  # 均方根值
    #     square_root_amplitude = (np.sum(np.sqrt(abs(window))) / cols) ** 2  # 方根幅值
    #     variance = np.var(window)  # 方差
    #     std = np.std(window)  # 标准差
    #     kurtosis = stats.kurtosis(window)[0]  # 峭度
    #     skewness = stats.skew(window)[0]  # 偏度
    # # mean_amplitude = np.sum(np.abs(data), axis=1) / cols  # 平均幅值 == 绝对平均值
    # # 无量纲统计量
    #     clearance_factor = peak_value / square_root_amplitude  # 裕度指标
    #     shape_factor = rms / abs_mean  # 波形指标
    #     impulse_factor = peak_value / abs_mean  # 脉冲指标
    #     crest_factor = peak_value / rms  # 峰值指标
    # # kurtosis_factor = kurtosis / (rms**4)  # 峭度指标
        features.append([max_value,peak_value, min_value, mean])

    return np.array(features, dtype=object)
#--------------------------------------频域特征提取------------------------------------------------
def frequencyDomainFeatures(src, fs ,window_size, step):
    # Calculate FFT of the signal
    pinyu = []
    for i in range(0, len(src) - window_size + 1, step):
        window = src[i:i + window_size]
        N = len(window)
        n = np.arange(N)
        freq = n * fs / N
        f = np.abs(np.fft.fft(window, N) * 2 / N)
        x = f[:N//2]  # Frequency amplitude values
        freq = freq[:N//2]  # Frequency values
    # Calculate frequency domain features
        frequencystructMF= np.mean(x)  # Mean frequency
        frequencystructFC = np.sum(freq.reshape(-1, 1) * x) / np.sum(x)# Center of mass frequency (centroid)
        #frequencystructRMSF = np.sqrt(np.sum(freq[:, np.newaxis] ** 2 * x) / np.sum(x))# Root mean square frequency
        #frequencystructRVF = np.sqrt(np.sum((freq[:, np.newaxis] - frequencystructFC)**2 * x) / np.sum(x))
        pinyu.append([frequencystructMF, frequencystructFC])
    return np.array(pinyu, dtype=object)
# frequencystructMF（Mean Frequency）：这是信号的平均频率。它计算了信号频谱中每个频率分量的加权平均值，其中权重是每个频率分量的幅度。这个特征可以用来表示信号的中心频率位置。
#
# frequencystructFC（Center of Mass Frequency）：也称为信号的质心频率或谱心频率。它是信号频谱的质心，表示在频域上信号的分布中心。计算方法是将每个频率分量的频率与其幅度相乘，然后将所有这些乘积相加，最后除以幅度的总和。这个特征可以用来确定信号的主要频率成分。
#
# frequencystructRMSF（Root Mean Square Frequency）：这是信号的频域上的均方根频率。它考虑了信号频谱中各频率分量的权重，频率的平方以及它们的幅度。计算方法是将每个频率分量的频率的平方与其幅度相乘，然后将所有这些乘积相加，最后取平方根。这个特征用于表示信号频谱的分布的宽度和离散性。
#
# frequencystructRVF（Root Variance Frequency）：这是信号频域上的根方差频率。它度量了信号频谱中各频率分量相对于信号质心频率的偏离程度。计算方法是将每个频率分量的频率与质心频率之差的平方与其幅度相乘，然后将所有这些乘积相加，最后取平方根。这个特征用于描述信号频谱的离散性和分布的扩散程度。
#
# 这些特征通常用于分析信号的频谱特性，以便更好地理解信号的频域信息和特征。
#--------------------------------------数据读取------------------------------------------------
def load_csv(path):
    data_read = pd.read_csv(path,index_col=False, encoding="utf-8", quoting=3,header = None)
    list = data_read.values.tolist()
    data = np.array(list)
    #print(data.shape)
    #print(data)
    return data
data = load_csv("个人数据集，本代码读取的是.CSV文件")
hangshu = data.shape[0]
lieshu  = data.shape[1]
print(hangshu,lieshu)
pitch =  np.array([data[:,0]])
pitch = pitch.T
roll =  np.array([data[:,1]])
roll = roll.T
yaw =  np.array([data[:,2]])
yaw = yaw.T
fs = 400
window_size = 4
step = 1
#step_size= 4
sampling_frequency=400
#提取时域特征
time1 = get_time_domain_feature(pitch, window_size, step)
pinyu1 = frequencyDomainFeatures(pitch,fs, window_size, step)
time2 = get_time_domain_feature(roll, window_size, step)
pinyu2 = frequencyDomainFeatures(roll,fs, window_size, step)
time3 = get_time_domain_feature(yaw, window_size, step)
pinyu3 = frequencyDomainFeatures(yaw,fs, window_size, step)
combined_time = [list(t1) + list(t2) + list(t3) for t1, t2, t3 in zip(time1, time2, time3)]
combined_fft = [list(p1) + list(p2) + list(p3) for p1, p2, p3 in zip(pinyu1, pinyu2, pinyu3)]
csv_file = 'shiyu.csv'
if os.path.exists('shiyu.csv'):
    os.remove('shiyu.csv')
# 将矩阵导出为CSV文件
with open(csv_file, 'w', newline='') as file:
    writer = csv.writer(file)
    writer.writerows(combined_time)
print(f"时域矩阵已成功导出为CSV文件：{csv_file}")
csv_file = 'piny.csv'
if os.path.exists('piny.csv'):
    os.remove('piny.csv')
# 将矩阵导出为CSV文件
with open(csv_file, 'w', newline='') as file:
    writer = csv.writer(file)
    writer.writerows(combined_fft)
print(f"频域矩阵已成功导出为CSV文件：{csv_file}")