diff --git a/exp1_basic.py b/exp1_basic.py index 49bf475..aaaf532 100644 --- a/exp1_basic.py +++ b/exp1_basic.py @@ -23,7 +23,7 @@ import reservoirpy.nodes as rpn rpy.verbosity(0) rpy.set_seed(42) -from tools import load_data +from tools import load_data, visualize_data_diff, visualize_psd_diff def build_esn_model(units, lr, sr, reg, output_dim, reservoir_node=None): @@ -50,28 +50,10 @@ def evaluate_results(model, clean_data, noisy_data, warmup=0, vis_data=True, vis # 3个维度分别绘制 noisy vs clean | prediction vs clean if vis_data: - plt.figure(figsize=(12, 6)) - for i in range(3): - plt.subplot(3, 2, 2*i+1) - plt.plot(clean_data[warmup:, i], label='clean') - plt.plot(noisy_data[warmup:, i], label='noisy') - plt.title(f"Dim {i+1} (Noisy vs Clean)") - plt.legend() - plt.subplot(3, 2, 2*i+2) - plt.plot(clean_data[warmup:, i], label='clean') - plt.plot(prediction[warmup:, i], label='prediction') - plt.title(f"Dim {i+1} (Prediction vs Clean)") - plt.legend() + visualize_data_diff(clean_data, noisy_data, prediction, warmup=0) # PSD 对比 if vis_psd: - plt.figure(figsize=(12, 6)) - for i in range(3): - plt.subplot(3, 1, i+1) - plt.psd(clean_data[warmup:, i], NFFT=1024, label='clean') - plt.psd(noisy_data[warmup:, i], NFFT=1024, label='noisy') - plt.psd(prediction[warmup:, i], NFFT=1024, label='prediction') - plt.title(f"Dim {i+1} PSD Comparison") - plt.legend() + visualize_psd_diff(clean_data, noisy_data, prediction, warmup=0) return snr, rmse, psnr, corr_coef, prediction diff --git a/exp2_comparison.py b/exp2_comparison.py index d1dd286..49508af 100644 --- a/exp2_comparison.py +++ b/exp2_comparison.py @@ -10,25 +10,203 @@ - 图表展示 (例如 Lyapunov 指数估计误差的折线图)。 """ -def load_data(): - # TODO: 生成并加载混沌系统数据（带噪和干净） - pass +import numpy as np +import matplotlib.pyplot as plt -def build_models(): - # TODO: 构建 RC 模型以及其他对比降噪模型（Wiener、卡尔曼、SVR、GBRT、LSTM等） - pass +from scipy.signal import wiener +from sklearn.ensemble import GradientBoostingRegressor -def train_models(): - # TODO: 训练各个模型，确保训练流程一致 - pass +import reservoirpy as rpy +import reservoirpy.nodes as rpn +rpy.verbosity(0) +rpy.set_seed(42) -def evaluate_models(): - # TODO: 评估不同模型的降噪效果和动力学指标差异 - pass +from tools import load_data, visualize_psd_diff, visualize_data_diff, print_exec_time + +def calculate_metrics(clean_data, prediction): + snr = 10 * np.log10(np.mean(clean_data**2) / np.mean((clean_data - prediction)**2)) + rmse = np.sqrt(np.mean((clean_data - prediction)**2)) + psnr = 20 * np.log10(np.max(clean_data) / rmse) + corr_coef = np.corrcoef(clean_data, prediction)[0, 1] + return snr, rmse, psnr, corr_coef + +@print_exec_time +def test_esn(clean_data, noisy_data, depth, **kwargs): + + # 分为训练集和测试集 按百分比 + train_size = int(len(clean_data) * 0.8) + train_clean_data = clean_data[:train_size] + train_noisy_data = noisy_data[:train_size] + test_clean_data = clean_data[train_size:] + test_noisy_data = noisy_data[train_size:] + # 形状均为 (n_samples, n_dimensions) -> (40000, 3) + + def build_esn_model(units, lr, sr, reg, output_dim, reservoir_node=None): + if reservoir_node is None: + reservoir_node = rpn.Reservoir(units=units, lr=lr, sr=sr) + output_node = rpn.Ridge(output_dim=output_dim, ridge=reg) + model = reservoir_node >> output_node + return model, reservoir_node + + def train_esn(model, clean_data, noisy_data, warmup=1000): + model.fit(noisy_data, clean_data, warmup=warmup) + prediction = model.run(noisy_data) + return model, prediction + + def validate_esn(model, noisy_data): + prediction = model.run(noisy_data) + return prediction + + results = [] + + lv1_units = kwargs.get('lv1_units', 300) + lv1_lr = kwargs.get('lv1_lr', 0.3) + lv1_sr = kwargs.get('lv1_sr', 0.9) + lv1_reg = kwargs.get('lv1_reg', 1e-6) + lv1_output_dim = kwargs.get('lv1_output_dim', 3) + lv1_model, reservoir_node = build_esn_model(lv1_units, lv1_lr, lv1_sr, lv1_reg, lv1_output_dim) + lv1_model, lv1_prediction = train_esn(lv1_model, train_clean_data, train_noisy_data) + lv1_validation = validate_esn(lv1_model, test_noisy_data) + lv1_metrics = calculate_metrics(test_clean_data, lv1_validation) + results.append(lv1_metrics) + + if depth > 1: + lvi_train_prediction = lv1_prediction + lvi_validation = lv1_validation + for i in range(2, depth+1): + lvi_units = kwargs.get(f'lv{i}_units', 300) + lvi_lr = kwargs.get(f'lv{i}_lr', 0.3) + lvi_sr = kwargs.get(f'lv{i}_sr', 0.9) + lvi_reg = kwargs.get(f'lv{i}_reg', 1e-6) + lvi_output_dim = kwargs.get(f'lv{i}_output_dim', 3) + lvi_model, _ = build_esn_model(lvi_units, lvi_lr, lvi_sr, lvi_reg, lvi_output_dim, reservoir_node) + lvi_model, lvi_train_prediction = train_esn(lvi_model, train_clean_data, lvi_train_prediction) + lvi_validation = validate_esn(lvi_model, lvi_validation) + lvi_metrics = calculate_metrics(test_clean_data, lvi_validation) + results.append(lvi_metrics) + + visualize_data_diff(test_clean_data, test_noisy_data, lvi_validation, warmup=0, title_prefix='ESN') + + return results, lvi_validation + +@print_exec_time +def test_kalman(clean_data, noisy_data, process_variance=1e-5, measurement_variance=1e-1): + # 假设 1D 卡尔曼模型: x(k+1) = x(k) + 过程噪声, z(k) = x(k) + 测量噪声 + # 对每个通道进行滤波 + # process_variance = 1e-5, 1e-4, 1e-3 + # measurement_variance = 1e-3, 1e-2, 1e-1 + # 1e-5, 1e-3 -> Kalman Metrics: (5.690448393281406 , 0.5213188005972011, 14.43781885284762 , 0.9999978828900199) + # 1e-3, 1e-1 -> Kalman Metrics: (5.6905031524143475, 0.5213155140166996, 14.437873611980564, 0.9999978828900199) + filtered_data = np.zeros_like(noisy_data) + for ch in range(noisy_data.shape[1]): + n_samples = noisy_data.shape[0] + x_est = 0.0 # 初始状态估计 + p_est = 1.0 # 初始估计协方差 + + for k in range(n_samples): + # 预测阶段 + x_pred = x_est + p_pred = p_est + process_variance + + # 更新阶段 + z_k = noisy_data[k, ch] # 测量值 + K = p_pred / (p_pred + measurement_variance) + x_est = x_pred + K * (z_k - x_pred) + p_est = (1 - K) * p_pred + + filtered_data[k, ch] = x_est + + metrics = calculate_metrics(clean_data, filtered_data) + return metrics, filtered_data + +@print_exec_time +def test_wiener(clean_data, noisy_data, mysize=35): + # mysize=[5, 7] + # mysize = 5 -> Wiener Metrics: (7.082912643590872 , 0.4440993882054787, 15.830283103157086, 0.9999978828900199) + # mysize = 7 -> Wiener Metrics: (8.393434517491048 , 0.3819038892573478, 17.14080497705726 , 0.9999978828900199) + # mysize = 9 -> Wiener Metrics: (9.396860794197918 , 0.3402379612434453, 18.144231253764133, 0.9999978828900199) + # mysize = 11 -> Wiener Metrics: (10.142598565197332, 0.3122452775883699, 18.889969024763545, 0.9999978828900199) + # mysize = 15 -> Wiener Metrics: (11.128377345656945, 0.2787449061849072, 19.87574780522316 , 0.9999978828900199) + # mysize = 25 -> Wiener Metrics: (11.145673899312722, 0.27819038279543434, 19.89304435887894, 0.9999978828900199) + filtered_data = np.zeros_like(noisy_data) + for i in range(noisy_data.shape[1]): + filtered_data[:, i] = wiener(noisy_data[:, i], mysize=mysize) + + metrics = calculate_metrics(clean_data, filtered_data) + return metrics, filtered_data + +@print_exec_time +def test_gbrt(clean_data, noisy_data, max_depth=13): + ''' + dep-13 -> GBRT Metrics: (5.187844922790017, 0.5523744230155584, 13.935215382356231, 0.9999978828900199) + dep-11 -> GBRT Metrics: (5.408985965603853, 0.5384885940824224, 14.156356425170069, 0.9999978828900199) + dep-9 -> GBRT Metrics: (5.298970192515195, 0.5453524862821368, 14.046340652081408, 0.9999978828900199) + dep-7 -> GBRT Metrics: (4.410332845010704, 0.604100457756635 , 13.15770330457692 , 0.9999978828900199) + ''' + gbrt = GradientBoostingRegressor(n_estimators=100, learning_rate=0.1, max_depth=max_depth, random_state=42) + filtered_data = np.zeros_like(noisy_data) # 初始化 filtered_data，形状与 noisy_data 相同 + + for i in range(noisy_data.shape[1]): # 遍历通道 + x_vals_all = np.arange(len(noisy_data)).reshape(-1, 1) # 完整数据集时间索引 + + gbrt.fit(x_vals_all, noisy_data[:, i]) # 在 *整个数据集* 上训练，目标是 noisy_data的 *单个通道* + filtered_data[:, i] = gbrt.predict(x_vals_all) # 预测 *整个数据集*，输入是 *时间索引* + + metrics = calculate_metrics(clean_data, filtered_data) + visualize_data_diff(clean_data, noisy_data, filtered_data, warmup=0, title_prefix='GBRT_FullData') # 修改 title_prefix 以区分 + return metrics, filtered_data + +@print_exec_time +def test_moving_average(clean_data, noisy_data, window_size=10): + # 10 -> Moving Average Metrics: (7.742861478461627, 0.4116069821396233, 16.49023193802784 , 0.9999978828900199) + # 15 -> Moving Average Metrics: (5.702762184458747, 0.5205802622393892, 14.450132644024961, 0.9999978828900199) + filtered_data = np.zeros_like(noisy_data) + + # 对每个通道应用移动平均滤波 + for i in range(noisy_data.shape[1]): + cumsum = np.cumsum(np.insert(noisy_data[:, i], 0, 0)) + avg = (cumsum[window_size:] - cumsum[:-window_size]) / window_size + # 从索引 window_size-1 开始赋值，确保形状一致 + filtered_data[window_size-1:, i] = avg + for j in range(window_size - 1): + filtered_data[j, i] = np.mean(noisy_data[0:j+1, i]) + + metrics = calculate_metrics(clean_data, filtered_data) + return metrics, filtered_data if __name__ == "__main__": - # 主对比实验流程 - load_data() - build_models() - train_models() - evaluate_models() + import argparse + parser = argparse.ArgumentParser() + parser.add_argument('--mode', type=str, default='esn', help='选择测试的降噪方法') + args = parser.parse_args() + + # 数据加载与预处理 + clean_data, noisy_data = load_data(system='lorenz', noise='gaussian', + intensity=0.8, init=[1, 1, 1], + n_timesteps=30000, transient=10000, h=0.01) + + if args.mode == 'gbrt' or args.mode == 'all': + gbrt_metrics, gbrt_prediction = test_gbrt(clean_data, noisy_data) + print("GBRT Metrics:", gbrt_metrics) + #visualize_data_diff(clean_data, noisy_data, gbrt_prediction, warmup=0, title_prefix='GBRT') + + if args.mode == 'kalman' or args.mode == 'all': + kalman_metrics, kalman_prediction = test_kalman(clean_data, noisy_data) + print("Kalman Metrics:", kalman_metrics) + visualize_data_diff(clean_data, noisy_data, kalman_prediction, warmup=0, title_prefix='Kalman') + + if args.mode == 'wiener' or args.mode == 'all': + wiener_metrics, wiener_prediction = test_wiener(clean_data, noisy_data) + print("Wiener Metrics:", wiener_metrics) + visualize_data_diff(clean_data, noisy_data, wiener_prediction, warmup=0, title_prefix='Wiener') + + if args.mode == 'esn' or args.mode == 'all': + esn_metrics, esn_prediction = test_esn(clean_data, noisy_data, depth=3) + print("ESN Metrics:", esn_metrics[-1]) + + if args.mode == 'moving_average' or args.mode == 'all': + ma_metrics, ma_prediction = test_moving_average(clean_data, noisy_data) + print("Moving Average Metrics:", ma_metrics) + visualize_data_diff(clean_data, noisy_data, ma_prediction, warmup=0, title_prefix='Moving Average') + + plt.show() \ No newline at end of file diff --git a/run_all_tests.bat b/run_all_tests.bat new file mode 100644 index 0000000..d092729 --- /dev/null +++ b/run_all_tests.bat @@ -0,0 +1,18 @@ +@echo off +rem 启动 ESN 测试 +start "" python "C:/Users/AresZz/Desktop/rc denoise/exp/exp2_comparison.py" --mode esn + +rem 启动 Kalman 测试 +start "" python "C:/Users/AresZz/Desktop/rc denoise/exp/exp2_comparison.py" --mode kalman + +rem 启动 Wiener 测试 +start "" python "C:/Users/AresZz/Desktop/rc denoise/exp/exp2_comparison.py" --mode wiener + +rem 启动 GBRT 测试 +start "" python "C:/Users/AresZz/Desktop/rc denoise/exp/exp2_comparison.py" --mode gbrt + +rem 启动移动平均测试 +start "" python "C:/Users/AresZz/Desktop/rc denoise/exp/exp2_comparison.py" --mode moving_average + +echo 测试已启动，等待所有测试完成后请按任意键退出... +pause diff --git a/tools.py b/tools.py index 0783e99..4b987b1 100644 --- a/tools.py +++ b/tools.py @@ -1,5 +1,15 @@ -from reservoirpy import datasets +import time import numpy as np +from matplotlib import pyplot as plt +from reservoirpy import datasets + +def print_exec_time(func): + def wrapper(*args, **kwargs): + start = time.time() + result = func(*args, **kwargs) + print(f"Function {func.__name__} executed in {time.time() - start:.4f} seconds.") + return result + return wrapper def add_noise(data, noise_type='gaussian', intensity=0.1, **kwargs): """ @@ -194,6 +204,32 @@ def load_data(system='lorenz', init='random', noise=None, intensity=0.1, h=0.01, return clean_data, noisy_data +def visualize_data_diff(clean_data, noisy_data, prediction, warmup=0, title_prefix=""): + plt.figure(figsize=(12, 6)) + for i in range(3): + plt.subplot(3, 2, 2*i+1) + plt.plot(noisy_data[warmup:, i], label='noisy') + plt.plot(clean_data[warmup:, i], label='clean') + plt.title(f"{title_prefix} Dim {i+1} (Noisy vs Clean)") + plt.legend() + plt.subplot(3, 2, 2*i+2) + plt.plot(prediction[warmup:, i], label='prediction') + plt.plot(clean_data[warmup:, i], label='clean') + plt.title(f"{title_prefix} Dim {i+1} (Prediction vs Clean)") + plt.legend() + plt.tight_layout() + +def visualize_psd_diff(clean_data, noisy_data, prediction, warmup=0, title_prefix=""): + plt.figure(figsize=(12, 6)) + for i in range(3): + plt.subplot(3, 1, i+1) + plt.psd(clean_data[warmup:, i], NFFT=1024, label='clean') + plt.psd(noisy_data[warmup:, i], NFFT=1024, label='noisy') + plt.psd(prediction[warmup:, i], NFFT=1024, label='prediction') + plt.title(f"{title_prefix} Dim {i+1} PSD Comparison") + plt.legend() + plt.tight_layout() + if __name__ == "__main__": # 测试数据加载和噪声添加 clean_data, noisy_data = load_data(system='lorenz', noise=['gaussian', 'colored'], intensity=0.1, n_timesteps=10000, transient=1000)