解决第一张图vmd各模态分解图为空白的问题import os
import json
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import IsolationForest
from scipy.signal import savgol_filter, find_peaks
from scipy.fft import fft, fftfreq
from sklearn.metrics import mean_squared_error
from vmdpy import VMD  # 需要安装vmdpy库: pip install vmdpy
import warnings

# 忽略不必要的警告
warnings.filterwarnings("ignore")

# 解决中文显示问题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False    # 解决负号显示问题

# 参数文件路径
PARAM_FILE = 'best_params.json'

def calculate_dominant_frequency(imf, sampling_rate):
    N = len(imf)
    yf = fft(imf)
    xf = fftfreq(N, d=1/sampling_rate)
    positive_freqs = xf[:N//2]
    positive_yf = np.abs(yf[:N//2])
    
    if len(positive_yf) == 0:
        return 0.0
    
    idx = np.argmax(positive_yf)
    dominant_freq = positive_freqs[idx]
    return dominant_freq

def vmd_spectral_decomposition(signal, time_seconds, alpha=2000, tau=0, K=5, DC=0, init=1):
    fs = 1 / np.mean(np.diff(time_seconds))  # 采样率
    imfs, u_hat, omega = VMD(signal, alpha, tau, K, DC, init, fs)

    num_IMFs = imfs.shape[0]
    print(f"分解得到的IMF数量: {num_IMFs}")

    sampling_rate = fs
    print(f"估计采样率: {sampling_rate:.4f} Hz")
    
    dominant_freqs = []
    for i, imf in enumerate(imfs):
        dom_freq = calculate_dominant_frequency(imf, sampling_rate)
        dominant_freqs.append(dom_freq)
        print(f"IMF {i+1}: 主频率 = {dom_freq:.4f} Hz")
    
    frequency_threshold = adaptive_frequency_threshold(dominant_freqs)
    print(f"自动确定的频率阈值为: {frequency_threshold:.4f} Hz")
    
    low_freq_IMF_indices = [i for i, freq in enumerate(dominant_freqs) if freq < frequency_threshold]
    print(f"需要去除的低频IMF编号: {low_freq_IMF_indices}")
    
    if low_freq_IMF_indices:
        removed_signal = np.sum(imfs[low_freq_IMF_indices, :], axis=0)
        print(f"Length of original signal: {len(signal)}")
        print(f"Length of removed signal: {len(removed_signal)}")
        
        # 确保 removed_signal 和 signal 的长度相同
        if len(signal) > len(removed_signal):
            removed_signal = np.pad(removed_signal, (0, len(signal) - len(removed_signal)), mode='edge')
        
        cleaned_signal = signal - removed_signal
    else:
        cleaned_signal = signal
    
    return cleaned_signal, low_freq_IMF_indices, imfs, frequency_threshold

def adaptive_frequency_threshold(dominant_freqs):
    hist, bin_edges = np.histogram(dominant_freqs, bins='auto', density=True)
    threshold_idx = np.argmax(hist)
    frequency_threshold = (bin_edges[threshold_idx] + bin_edges[threshold_idx + 1]) / 2
    return frequency_threshold

def dynamic_window_size(cleaned_signal):
    autocorrelation = np.correlate(cleaned_signal, cleaned_signal, mode='full')
    autocorrelation = autocorrelation[autocorrelation.size // 2:]
    
    half_max = autocorrelation[0] * 0.5
    peaks, _ = find_peaks(-autocorrelation, height=-half_max)
    if len(peaks) > 0:
        window_size = peaks[0]
    else:
        window_size = min(24, len(cleaned_signal) // 2)
    return max(12, window_size)

def moving_average_forecast(series, window_size):
    forecast = []
    for t in range(len(series) - window_size + 1):
        forecast.append(series[t:t + window_size].mean())
    return np.array(forecast)

def plot_imfs(time, IMFs, low_freq_indices):
    num_IMFs = IMFs.shape[0]
    plt.figure(figsize=(12, 3*num_IMFs))
    
    if len(time) != IMFs.shape[1]:
        print("Warning: length of time and IMFs do not match.")
        return
    
    for i in range(num_IMFs):
        plt.subplot(num_IMFs, 1, i + 1)
        plt.plot(time, IMFs[i], label=f'IMF {i+1}', color='b' if i not in low_freq_indices else 'r')
        if i in low_freq_indices:
            plt.legend(loc='upper right')
            plt.title(f'IMF {i+1} (低频, 温度影响)')
        else:
            plt.legend(loc='upper right')
            plt.title(f'IMF {i+1}')
        plt.xlabel('时间')
        plt.ylabel('挠度')
    
    plt.tight_layout()
    plt.show()

def save_best_params(params):
    with open(PARAM_FILE, 'w', encoding='utf-8') as f:
        json.dump(params, f, ensure_ascii=False, indent=4)
    print(f"最佳参数已保存到 {PARAM_FILE}。")

def load_best_params():
    if not os.path.exists(PARAM_FILE):
        print(f"参数文件 {PARAM_FILE} 不存在。")
        return None
    with open(PARAM_FILE, 'r', encoding='utf-8') as f:
        params = json.load(f)
    print(f"已加载最佳参数: {params}")
    return params

def estimate_contamination(series):
    mean = np.mean(series)
    std = np.std(series)
    threshold_upper = mean + 3 * std
    threshold_lower = mean - 3 * std
    num_anomalies = np.sum((series > threshold_upper) | (series < threshold_lower))
    contamination = num_anomalies / len(series)
    contamination = min(max(contamination, 0.01), 0.05)
    print(f"自动估计的异常比例（contamination）: {contamination:.4f}")
    return contamination

def main_automated_process(file_path):
    try:
        df = pd.read_excel(file_path)
    except Exception as e:
        print(f"读取Excel文件时出错: {e}")
        return None, None
    
    print("Excel文件读取成功。")
    print(f"数据总行数: {len(df)}")
    print(f"列名: {df.columns.tolist()}")
    
    required_columns = ['时间', '挠度']
    for col in required_columns:
        if col not in df.columns:
            print(f"Excel文件中缺少 '{col}' 列，请检查列名。")
            return None, None
    
    try:
        df['时间'] = pd.to_datetime(df['时间'])
    except Exception as e:
        print(f"转换 '时间' 列为 datetime 格式时出错: {e}")
        return None, None
    
    window_size = 11
    poly_order = 3
    if window_size % 2 == 0:
        window_size += 1
        print(f"window_size为偶数，已调整为 {window_size}")
    
    try:
        df['平滑挠度'] = savgol_filter(df['挠度'], window_size, poly_order)
    except Exception as e:
        print(f"Savgol 滤波器应用时出错: {e}")
        return None, None
    
    contamination = estimate_contamination(df['平滑挠度'].values)
    
    iso_forest = IsolationForest(contamination=contamination, random_state=42)
    try:
        anomalies = iso_forest.fit_predict(df[['平滑挠度']])
    except Exception as e:
        print(f"孤立森林异常检测时出错: {e}")
        return None, None
    
    df['挠度_去除异常'] = df['平滑挠度']
    df.loc[anomalies == -1, '挠度_去除异常'] = np.nan
    num_anomalies = (anomalies == -1).sum()
    print(f"检测到的异常值数量: {num_anomalies}")
    
    cleaned_data = df[['时间', '挠度_去除异常']].dropna()
    
    if cleaned_data.empty:
        print("所有数据点都被标记为异常，无法继续处理。")
        return df, None
    
    try:
        time_seconds = (cleaned_data['时间'] - cleaned_data['时间'].iloc[0]).dt.total_seconds().values
    except Exception as e:
        print(f"转换时间列为秒数时出错: {e}")
        return None, None
    
    try:
        cleaned_signal, removed_IMFs, IMFs, frequency_threshold = vmd_spectral_decomposition(
            signal=cleaned_data['挠度_去除异常'].values,
            time_seconds=time_seconds
        )
    except Exception as e:
        print(f"VMD分解过程中出错: {e}")
        return df, None

    df.loc[cleaned_data.index[:len(cleaned_signal)], '去除温致挠度'] = cleaned_signal

    plot_imfs(cleaned_data['时间'].iloc[:len(cleaned_signal)], IMFs, removed_IMFs)
    
    window_size_forecast = dynamic_window_size(cleaned_signal)
    print(f"动态确定的移动平均窗口大小 (window_size_forecast): {window_size_forecast}")
    
    try:
        forecast = moving_average_forecast(cleaned_signal, window_size_forecast)
        if len(forecast) == 0:
            print("移动平均预测生成的预测数据为空。")
        else:
            print(f"Length of cleaned_signal: {len(cleaned_signal)}")
            print(f"Length of forecast: {len(forecast)}")
            actual = cleaned_signal[window_size_forecast - 1:]
            predicted = forecast[:len(actual)]
            rmse = np.sqrt(mean_squared_error(actual, predicted))
            print(f"移动平均预测的RMSE: {rmse:.4f}")
            forecast_series = np.concatenate([cleaned_signal[:window_size_forecast - 1], forecast])
            df.loc[cleaned_data.index[:len(forecast_series)], '预测'] = forecast_series
    except Exception as e:
        print(f"移动平均预测过程中出错: {e}")
    
    used_params = {
        'contamination': contamination,
        'frequency_threshold': frequency_threshold,
        'window_size_forecast': window_size_forecast
    }
    save_best_params(used_params)
    
    try:
        plt.figure(figsize=(14, 8))
        
        plt.subplot(2, 1, 1)
        plt.plot(cleaned_data['时间'].iloc[:len(cleaned_signal)], cleaned_data['挠度_去除异常'].iloc[:len(cleaned_signal)], label='原始信号', color='b', alpha=0.5)
        plt.plot(cleaned_data['时间'].iloc[:len(cleaned_signal)], cleaned_signal, label='去除低频IMF后的信号', color='g')
        plt.legend(loc='upper right')
        plt.title('原始信号 vs 去除低频IMF后的信号')
        plt.xlabel('时间')
        plt.ylabel('挠度')
        
        if '预测' in df.columns:
            plt.subplot(2, 1, 2)
            plt.plot(cleaned_data['时间'].iloc[:len(forecast_series)], cleaned_signal[:len(forecast_series)], label='去除低频IMF后的信号', color='g')
            plt.plot(cleaned_data['时间'].iloc[:len(forecast_series)], forecast_series, label='移动平均预测', color='r')
            plt.legend(loc='upper right')
            plt.title('去除低频IMF后的信号及其预测')
            plt.xlabel('时间')
            plt.ylabel('挠度')
        plt.tight_layout()
        plt.show()
    except Exception as e:
        print(f"可视化过程中出错: {e}")
    
    return df, used_params

def load_and_process_with_best_params(file_path, best_params):
    contamination = best_params.get('contamination', 0.02)
    frequency_threshold = best_params.get('frequency_threshold', 0.1)
    window_size_forecast = best_params.get('window_size_forecast', 24)
    
    try:
        df = pd.read_excel(file_path)
    except Exception as e:
        print(f"读取Excel文件时出错: {e}")
        return None
    
    print("Excel文件读取成功。")
    print(f"数据总行数: {len(df)}")
    print(f"列名: {df.columns.tolist()}")
    
    required_columns = ['时间', '挠度']
    for col in required_columns:
        if col not in df.columns:
            print(f"Excel文件中缺少 '{col}' 列，请检查列名。")
            return None
    
    try:
        df['时间'] = pd.to_datetime(df['时间'])
    except Exception as e:
        print(f"转换 '时间' 列为 datetime 格式时出错: {e}")
        return None
    
    window_size = 11
    poly_order = 3
    if window_size % 2 == 0:
        window_size += 1
        print(f"window_size为偶数，已调整为 {window_size}")
    
    try:
        df['平滑挠度'] = savgol_filter(df['挠度'], window_size, poly_order)
    except Exception as e:
        print(f"Savgol 滤波器应用时出错: {e}")
        return None
    
    iso_forest = IsolationForest(contamination=contamination, random_state=42)
    try:
        anomalies = iso_forest.fit_predict(df[['平滑挠度']])
    except Exception as e:
        print(f"孤立森林异常检测时出错: {e}")
        return None
    
    df['挠度_去除异常'] = df['平滑挠度']
    df.loc[anomalies == -1, '挠度_去除异常'] = np.nan
    num_anomalies = (anomalies == -1).sum()
    print(f"检测到的异常值数量: {num_anomalies}")
    
    cleaned_data = df[['时间', '挠度_去除异常']].dropna()
    
    if cleaned_data.empty:
        print("所有数据点都被标记为异常，无法继续处理。")
        return df
    
    try:
        time_seconds = (cleaned_data['时间'] - cleaned_data['时间'].iloc[0]).dt.total_seconds().values
    except Exception as e:
        print(f"转换时间列为秒数时出错: {e}")
        return None
    
    try:
        cleaned_signal, removed_IMFs, IMFs, _ = vmd_spectral_decomposition(
            signal=cleaned_data['挠度_去除异常'].values,
            time_seconds=time_seconds
        )
    except Exception as e:
        print(f"VMD分解过程中出错: {e}")
        return df
    
    df.loc[cleaned_data.index[:len(cleaned_signal)], '去除温致挠度'] = cleaned_signal

    plot_imfs(cleaned_data['时间'].iloc[:len(cleaned_signal)], IMFs, removed_IMFs)
    
    try:
        forecast = moving_average_forecast(cleaned_signal, window_size_forecast)
        if len(forecast) == 0:
            print("移动平均预测生成的预测数据为空。")
        else:
            print(f"Length of cleaned_signal: {len(cleaned_signal)}")
            print(f"Length of forecast: {len(forecast)}")
            actual = cleaned_signal[window_size_forecast - 1:]
            predicted = forecast[:len(actual)]
            rmse = np.sqrt(mean_squared_error(actual, predicted))
            print(f"移动平均预测的RMSE: {rmse:.4f}")
            forecast_series = np.concatenate([cleaned_signal[:window_size_forecast - 1], forecast])
            df.loc[cleaned_data.index[:len(forecast_series)], '预测'] = forecast_series
    except Exception as e:
        print(f"移动平均预测过程中出错: {e}")
    
    try:
        plt.figure(figsize=(14, 8))
        
        plt.subplot(2, 1, 1)
        plt.plot(cleaned_data['时间'].iloc[:len(cleaned_signal)], cleaned_data['挠度_去除异常'].iloc[:len(cleaned_signal)], label='原始信号', color='b', alpha=0.5)
        plt.plot(cleaned_data['时间'].iloc[:len(cleaned_signal)], cleaned_signal, label='去除低频IMF后的信号', color='g')
        plt.legend(loc='upper right')
        plt.title('原始信号 vs 去除低频IMF后的信号')
        plt.xlabel('时间')
        plt.ylabel('挠度')
        
        if '预测' in df.columns:
            plt.subplot(2, 1, 2)
            plt.plot(cleaned_data['时间'].iloc[:len(forecast_series)], cleaned_signal[:len(forecast_series)], label='去除低频IMF后的信号', color='g')
            plt.plot(cleaned_data['时间'].iloc[:len(forecast_series)], forecast_series, label='移动平均预测', color='r')
            plt.legend(loc='upper right')
            plt.title('去除低频IMF后的信号及其预测')
            plt.xlabel('时间')
            plt.ylabel('挠度')
        plt.tight_layout()
        plt.show()
    except Exception as e:
        print(f"可视化过程中出错: {e}")
    
    return df, best_params

def display_menu():
    print("\n请选择操作模式:")
    print("1. 自适应选择参数 (自动参数优化)")
    print("2. 使用上次的参数")
    print("3. 退出")
    choice = input("请输入选项 (1/2/3): ").strip()
    return choice

def main():
    file_path = './挠度数据910.xlsx'  # 替换为实际文件路径
    
    while True:
        choice = display_menu()
        
        if choice == '1':
            print("\n您选择了自适应选择参数 (自动参数优化)。")
            df_result, optimal_params = main_automated_process(file_path)
            if optimal_params:
                print("\n自动化处理完成。")
                print(f"使用的参数组合: {optimal_params}")
        elif choice == '2':
            print("\n您选择了使用上次的参数。")
            best_params = load_best_params()
            if best_params is None:
                print("未找到上次的参数，请先进行自适应选择参数 (选项1)。")
                continue
            df_result, _ = load_and_process_with_best_params(file_path, best_params)
            if df_result is not None:
                print("\n使用上次的参数处理完成。")
        elif choice == '3':
            print("退出程序。")
            break
        else:
            print("无效的选项，请重新选择。")

if __name__ == "__main__":
    main()您选择了自适应选择参数 (自动参数优化)。
Excel文件读取成功。
数据总行数: 83557
列名: ['时间', '挠度', '温度']
自动估计的异常比例（contamination）: 0.0100
检测到的异常值数量: 822
分解得到的IMF数量: 5
估计采样率: 0.9576 Hz
IMF 1: 主频率 = 0.0000 Hz
IMF 2: 主频率 = 0.0339 Hz
IMF 3: 主频率 = 0.0526 Hz
IMF 4: 主频率 = 0.2846 Hz
IMF 5: 主频率 = 0.3790 Hz
自动确定的频率阈值为: 0.0474 Hz
需要去除的低频IMF编号: [0, 1]
Length of original signal: 82735
Length of removed signal: 82734
Warning: length of time and IMFs do not match.
动态确定的移动平均窗口大小 (window_size_forecast): 12
Length of cleaned_signal: 82735
Length of forecast: 82724
移动平均预测的RMSE: 1.3473
最佳参数已保存到 best_params.json。

自动化处理完成。
使用的参数组合: {'contamination': 0.01, 'frequency_threshold': 0.0473755787037037, 'window_size_forecast': 12}