``` def calculate_snr(signal, noise): """Calculate the SNR between a signal and its associated noise.""" if len(signal) != len(noise): raise ValueError("Signal and noise must have the same length.") signal_power = np.mean(np.abs(signal) ** 2) noise_power = np.mean(np.abs(noise) ** 2) if noise_power == 0: return float('inf') return 10 * np.log10(signal_power / noise_power) def add_noise_with_snr(original_signal, target_snr_db): """ Add Gaussian white noise to an input signal so that it reaches a specified SNR. :param original_signal: Original clean signal array. :param target_snr_db: Target SNR value in decibels. :return: Tuple containing noisy signal and added noise arrays. """ signal_power = np.mean(np.abs(original_signal) ** 2) snr_linear = 10 ** (target_snr_db / 10.0) noise_power = signal_power / snr_linear noise_std = np.sqrt(noise_power) noise = np.random.normal(0, noise_std, size=len(original_signal)) noisy_signal = original_signal + noise return noisy_signal, noise def wavelet_denoising(signal, level=3): """Perform Wavelet Transform based denoising on a given signal.""" wp = pywt.WaveletPacket(data=signal, wavelet='db1', mode='symmetric', maxlevel=level) nodes = [node.path for node in wp.get_level(level, 'natural')] thresholded_coeffs = {} for n in nodes: coeffs = wp[n].data sigma = np.median(np.abs(coeffs)) / 0.6745 threshold = sigma * np.sqrt(2 * np.log(len(coeffs))) coeffs[np.abs(coeffs) < threshold] = 0 # Soft-thresholding thresholded_coeffs[n] = coeffs.flatten() reconstructions = [] for k, v in sorted(thresholded_coeffs.items()): packet = pywt.WaveletPacket(None, wavelet='db1', mode='symmetric') packet[k] = v reconstructed = packet.reconstruct(update=False).real reconstructions.append(reconstructed) final_reconstruction = sum(reconstructions) / len(nodes) return final_reconstruction[:len(signal)] # Ensure output matches input length if __name__ == "__main__": # Generate test signals t = np.linspace(0, 1, num=2 ** 8) f1, f2 = 7, 13 # Frequencies of sinusoidal components clean_signal = np.sin(2 * np.pi * f1 * t) + np.cos(2 * np.pi * f2 * t) target_snrs = [5, 10, 20] results = {} fig, axs = plt.subplots(nrows=len(target_snrs), ncols=1, figsize=(10, 6 * len(target_snrs)), sharex=True) # 确保axs总是可以按索引访问的形式 if not isinstance(axs, np.ndarray): axs = [axs] for idx, target_snr_db in enumerate(target_snrs): noisy_signal, _ = add_noise_with_snr(clean_signal, target_snr_db) denoised_signal = wavelet_denoising(noisy_signal, level=3) actual_snr_before = calculate_snr(clean_signal, noisy_signal - clean_signal) actual_snr_after = calculate_snr(clean_signal, denoised_signal - clean_signal) print(f"\nResults for Target SNR={target_snr_db:.0f} dB:") print(f"- Actual Noisy SNR: {actual_snr_before:.2f} dB") print(f"- Actual Denoised SNR: {actual_snr_after:.2f} dB\n") ax = axs[idx] ax.plot(t, clean_signal, label="Original Signal", color='black', alpha=0.9) ax.plot(t, noisy_signal, label=f"Noisy Signal ({actual_snr_before:.2f} dB)", ls="--", lw=0.8) ax.plot(t, denoised_signal, label=f"Denoised Signal ({actual_snr_after:.2f} dB)", color='red') ax.set_title(f'Target SNR = {target_snr_db} dB') ax.grid(True) ax.legend(loc='upper right') plt.tight_layout(pad=3.0) plt.xlabel('Time (seconds)') plt.show()```如何用西交数据集验证去噪效果

<think>好的，我现在需要帮助用户理解如何用西交数据集验证他们提供的去噪代码的效果。首先，我要仔细阅读用户提供的Python代码，弄清楚各个函数的作用和整个流程。然后，我需要考虑如何将现有的测试信号替换为西交数据集的数据，并进行相应的调整。 首先，用户给出的代码包括三个主要函数：计算信噪比（SNR）、添加指定SNR的高斯白噪声、以及基于小波变换的去噪方法。主程序部分生成测试信号，添加不同SNR的噪声，去噪后计算并比较SNR，同时绘制信号对比图。用户现在想用西交数据集来验证去噪效果，所以需要将现有的测试信号替换为西交的数据。 接下来，我需要了解西交数据集的结构和格式。通常，这类数据集可能包含不同工况下的振动信号，可能有CSV或MAT文件，包含正常和故障状态的数据。假设用户的数据集是类似的，我需要指导他们如何加载这些数据，替换代码中的clean_signal。 首先，数据加载部分可能需要使用pandas读取CSV或scipy读取MAT文件。例如，用pd.read_csv加载数据，并提取所需的信号列。用户需要注意数据的长度和采样率，可能需要调整时间轴t的生成，以匹配实际信号的时长和采样点数量。 然后，替换原有的clean_signal生成代码。原来的clean_signal是合成的正弦波组合，现在需要替换为从数据集中加载的真实信号。用户需要确保加载的信号是一维数组，与噪声添加和去噪函数兼容。 接下来，可能需要对数据进行预处理。例如，数据可能需要进行归一化，因为不同的传感器数据可能有不同的幅值范围，归一化后可以更公平地比较去噪效果。此外，数据可能存在基线漂移或其他噪声，可能需要先进行初步的预处理，但用户可能希望保留这些真实噪声，以测试去噪算法的鲁棒性。 在添加噪声的部分，用户的add_noise_with_snr函数已经可以处理任意输入信号，所以替换数据后这部分应该不需要修改。但需要注意，如果数据集的采样率与原有代码中的t变量相关，可能需要调整时间轴，尤其是在绘图时正确显示时间单位。 去噪函数wavelet_denoising使用的是小波包变换，可能需要根据实际数据的特性调整小波基函数和分解层数。例如，db1小波适用于某些类型的信号，但对于更复杂的数据，可能需要尝试db4或其他小波，这需要用户根据实际情况调整参数。 验证部分，用户需要比较去噪前后的SNR。但这里有个问题：西交数据集中的信号本身可能已经包含噪声，因此要获取“干净”的信号作为基准可能困难。如果数据集包含正常状态下的信号作为干净信号，而故障信号作为噪声，或者有添加已知噪声的情况，那么可以计算SNR。否则，可能需要采用其他评估指标，如信噪比改善量（ΔSNR），或者使用视觉对比和时频分析来评估去噪效果。 此外，用户可能需要考虑使用不同的评估指标，如均方误差（MSE）、峰值信噪比（PSNR）、结构相似性指数（SSIM）等，特别是在没有干净基准的情况下。不过根据现有代码，用户可能希望继续使用SNR，因此需要确保数据集中有对应的干净和噪声信号。 可能还需要调整绘图部分，确保显示的数据量适中，避免过多的数据点导致图像不清晰。另外，如果数据较长，可能需要分段处理，或者选择代表性的片段进行可视化。 最后，用户需要运行修改后的代码，观察去噪后的SNR提升情况，并通过图表对比原始信号、噪声信号和去噪信号。如果效果不理想，可能需要调整小波去噪的参数，如小波类型、分解层数、阈值方法等，或者尝试其他去噪算法进行比较。 总结步骤应该是：加载数据→预处理→添加噪声（可选）→去噪→计算SNR→可视化。如果数据本身含有真实噪声，则无需添加，直接去噪并评估。但根据用户提供的代码，他们可能希望控制添加的噪声水平，因此可能需要将数据集中的干净信号作为基准，添加已知噪声，再测试去噪效果。 可能遇到的问题包括数据格式不匹配、信号长度不一致、采样率不同导致的时轴错误，以及小波参数不适合实际数据特征。用户需要逐一检查这些方面，确保代码适配数据集的特点。</think>以下是将代码适配到西交轴承数据集的实现步骤和注意事项，采用分步说明： $$ \text{SNR}_{\text{dB}} = 10 \log_{10}\left(\frac{P_{\text{signal}}}{P_{\text{noise}}}\right) $$ 一、数据准备阶段 1. 数据加载 ```python import pandas as pd # 假设数据集为CSV格式，包含振动信号列 def load_xjtu_data(file_path): df = pd.read_csv(file_path) # 假设'vibration'列为振动信号，'condition'为工况标签 raw_signal = df['vibration'].values # 转换为numpy数组 return raw_signal, df['condition'].iloc[0] # 返回信号和工况标签 ``` 2. 信号预处理 ```python def preprocess_signal(signal): # 标准化处理（保持与原始代码兼容性） signal = (signal - np.mean(signal)) / np.std(signal) # 低通滤波（可选，根据实际需求） from scipy.signal import butter, filtfilt b, a = butter(4, 0.1, 'lowpass') return filtfilt(b, a, signal) ``` 二、核心代码修改 ```python if __name__ == "__main__": # 加载西交数据集 xjtu_signal, condition = load_xjtu_data("Bearing1_1.csv") # 示例文件路径 clean_signal = preprocess_signal(xjtu_signal) # 时间轴生成（根据实际采样频率调整） sampling_rate = 25600 # 假设采样率25.6kHz t = np.arange(len(clean_signal)) / sampling_rate # 实际时间轴 ``` 三、评估指标增强 1. 增加时频分析 ```python from scipy.signal import spectrogram def plot_spectrogram(ax, signal, fs): f, t_spec, Sxx = spectrogram(signal, fs=fs, nperseg=256) ax.pcolormesh(t_spec, f, 10 * np.log10(Sxx), shading='gouraud') ax.set_ylabel('Frequency [Hz]') ``` 2. 增加多维度评估 ```python def calculate_rmse(original, denoised): return np.sqrt(np.mean((original - denoised)**2)) ``` 四、可视化优化 ```python # 在绘图部分增加时频分析 fig = plt.figure(figsize=(16, 12)) gs = gridspec.GridSpec(2, len(target_snrs), height_ratios=[3,1]) for idx, target_snr_db in enumerate(target_snrs): # 时域图 ax_time = plt.subplot(gs[0, idx]) # 原有时域绘图代码 # 频域图 ax_freq = plt.subplot(gs[1, idx]) plot_spectrogram(ax_freq, denoised_signal, sampling_rate) ``` 五、参数优化建议 1. 小波基选择 ```python # 尝试不同小波基函数 wavelet_types = ['db4', 'sym5', 'coif3'] # 机械信号常用小波基 for wavelet in wavelet_types: wp = pywt.WaveletPacket(data=signal, wavelet=wavelet, ...) ``` 2. 自适应阈值改进 ```python # 改进阈值计算方法 def adaptive_threshold(coeffs): N = len(coeffs) sigma = np.median(np.abs(coeffs)) / 0.6745 return sigma * np.sqrt(2 * np.log(N)) * (1 + 0.2 * np.log(N/1024)) ``` 六、完整验证流程 1. 数据加载 → 2. 预处理 → 3. 加噪 → 4. 去噪 → 5. 评估 → 6. 可视化 注意事项： 1. 采样率一致性：确保生成的时间轴与实际采样率匹配 2. 工况对应：不同故障类型需要分别验证（内圈/外圈/滚动体故障） 3. 数据长度：建议截取2^n长度数据（如4096点）以提高小波计算效率 4. 实时性验证：对于长信号可采用滑动窗口处理 典型输出结果解读： - SNR提升3-8dB为有效去噪 - 时频图中高频噪声成分应明显减少 - 故障特征频率（如BPFO/BPPI）在去噪后应更清晰 扩展建议： 1. 添加实际工况噪声：将高斯白噪声替换为实际轴承噪声 2. 迁移学习：使用西交数据训练去噪参数，应用到其他工业场景 3. 结合包络分析：对去噪后的信号进行包络谱分析以增强故障特征