当前位置：首页 > news >正文

用Python复现MATLAB经典案例：手把手教你处理温度传感器数据与消除60Hz工频干扰

news 2026/6/1 11:38:01

Python实战：温度传感器数据处理与60Hz工频干扰消除全指南

当波士顿机场的温度传感器每小时产生一组读数，或是电路中的开环电压信号被60Hz电源噪声污染时，传统MATLAB用户的第一反应可能是打开Signal Processing Toolbox。但在这个Python主导的开源时代，如何用NumPy、SciPy和Pandas这些工具实现同等甚至更优的信号处理效果？本文将彻底拆解两个经典案例，带您完成从MATLAB思维到Python实践的完美迁移。

1. 环境准备与数据加载

在开始信号处理前，我们需要配置合适的Python环境。推荐使用Anaconda创建独立环境：

conda create -n signal_processing python=3.9 conda activate signal_processing pip install numpy scipy pandas matplotlib seaborn

对于温度传感器数据，我们使用Pandas进行专业化的时间序列处理。假设原始数据来自CSV文件，包含时间戳和温度两列：

import pandas as pd import numpy as np from scipy import signal import matplotlib.pyplot as plt # 模拟波士顿机场温度数据（2011年1月） date_range = pd.date_range('2011-01-01', '2011-01-31 23:00:00', freq='H') np.random.seed(42) base_temp = 10 * np.sin(2*np.pi * date_range.dayofyear/365) # 季节趋势 daily_cycle = 5 * np.sin(2*np.pi * date_range.hour/24) # 昼夜周期 noise = np.random.normal(0, 1, len(date_range)) # 随机噪声 tempC = base_temp + daily_cycle + noise temp_data = pd.DataFrame({ 'datetime': date_range, 'tempC': tempC }).set_index('datetime')

提示：实际工程中建议使用pd.read_csv()加载真实传感器数据，并检查缺失值。使用df.interpolate()进行合理插值补全。

2. 温度数据的趋势提取与周期分析

2.1 移动平均滤波实现

24小时移动平均是提取日趋势的基础方法。与MATLAB的filter函数不同，Python中我们有多种实现选择：

# 方法1：使用NumPy的convolve window_size = 24 coeff_24h = np.ones(window_size)/window_size avg_24h_temp = np.convolve(tempC, coeff_24h, mode='same') # 方法2：使用Pandas的rolling（更适合时间序列） temp_data['24h_avg'] = temp_data['tempC'].rolling( window=24, center=True, min_periods=12).mean() # 可视化对比 plt.figure(figsize=(12,6)) plt.plot(temp_data.index, temp_data['tempC'], alpha=0.5, label='原始数据') plt.plot(temp_data.index, avg_24h_temp, 'r', label='NumPy移动平均') plt.plot(temp_data.index, temp_data['24h_avg'], 'g--', label='Pandas移动平均') plt.legend(); plt.ylabel('温度(℃)'); plt.title('24小时移动平均效果对比')

关键参数说明：

mode='same'：保证输出长度与输入一致
center=True：使窗口居中对齐，避免相位偏移
min_periods：设置最小计算窗口，处理边界情况

2.2 昼夜温差特征提取

通过移除日趋势，我们可以聚焦分析昼夜温差模式：

# 计算温差并重塑为(天数, 24小时)的矩阵 delta_temp = temp_data['tempC'] - temp_data['24h_avg'] delta_matrix = delta_temp.values.reshape(-1, 24) # 计算各小时的平均温差 hourly_effect = pd.DataFrame({ 'hour': range(24), 'mean_diff': delta_matrix.mean(axis=0), 'std_dev': delta_matrix.std(axis=0) }) # 可视化昼夜温差模式 plt.figure(figsize=(10,5)) plt.bar(hourly_effect['hour'], hourly_effect['mean_diff'], yerr=hourly_effect['std_dev']) plt.xlabel('一天中的小时'); plt.ylabel('与日均温差(℃)') plt.title('波士顿机场1月份昼夜温差特征')

3. 高级滤波技术实战

3.1 Savitzky-Golay滤波器应用

对于需要保留信号特征的平滑处理，SG滤波器是理想选择。Python中通过scipy.signal.savgol_filter实现：

# 不同阶数SG滤波器对比 window_length = 7 # 必须为奇数 temp_data['SG_cubic'] = signal.savgol_filter(tempC, window_length, 3) temp_data['SG_quartic'] = signal.savgol_filter(tempC, window_length, 4) # 效果可视化 plt.figure(figsize=(12,6)) plt.plot(temp_data.index, temp_data['tempC'], alpha=0.3, label='原始数据') plt.plot(temp_data.index, temp_data['SG_cubic'], 'r', label='三次SG滤波') plt.plot(temp_data.index, temp_data['SG_quartic'], 'g', label='四次SG滤波') plt.xlim(['2011-01-03','2011-01-05']) # 聚焦三天数据 plt.legend(); plt.title('Savitzky-Golay滤波器效果对比')

SG滤波器参数选择建议：

窗口长度：通常5-25之间的奇数，取决于信号特征
多项式阶数：一般2-4阶，过高会导致过拟合

3.2 中值滤波与异常值处理

当数据存在脉冲噪声时，中值滤波比均值滤波更具鲁棒性：

# 模拟含异常值的温度数据 temp_with_outliers = tempC.copy() outlier_indices = np.random.choice(len(tempC), size=20, replace=False) temp_with_outliers[outlier_indices] += np.random.uniform(-10,10,20) # 应用中值滤波 median_filtered = signal.medfilt(temp_with_outliers, kernel_size=5) # 可视化对比 plt.figure(figsize=(12,6)) plt.plot(temp_data.index, temp_with_outliers, 'o', markersize=3, label='含异常值') plt.plot(temp_data.index, median_filtered, 'r', linewidth=2, label='中值滤波') plt.legend(); plt.title('中值滤波对异常值的处理效果')

4. 60Hz工频干扰消除方案

4.1 电源噪声的模拟与采集

我们先模拟一个受60Hz干扰的开环电压信号：

fs = 1000 # 采样率1kHz t = np.arange(0, 1, 1/fs) signal_freq = 10 # 10Hz有用信号 true_voltage = 2 * np.sin(2*np.pi * signal_freq * t) noise_60hz = 0.5 * np.sin(2*np.pi * 60 * t + np.pi/3) measured_voltage = true_voltage + noise_60hz + np.random.normal(0,0.1,fs) plt.figure(figsize=(12,4)) plt.plot(t, measured_voltage, label='实测信号') plt.plot(t, true_voltage, 'r', linewidth=2, label='真实信号') plt.xlabel('时间(s)'); plt.ylabel('电压(V)') plt.legend(); plt.title('受60Hz干扰的电压信号')

4.2 自适应陷波滤波器设计

针对固定频率干扰，陷波滤波器是最直接有效的解决方案：

# 设计50Hz陷波滤波器 def notch_filter(fs, f0, Q=30): nyq = 0.5 * fs w0 = f0 / nyq b, a = signal.iirnotch(w0, Q) return b, a b, a = notch_filter(fs=1000, f0=60, Q=30) filtered_voltage = signal.filtfilt(b, a, measured_voltage) # 频域分析 f, Pxx = signal.welch(measured_voltage, fs, nperseg=1024) f_filt, Pxx_filt = signal.welch(filtered_voltage, fs, nperseg=1024) plt.figure(figsize=(12,8)) plt.subplot(2,1,1) plt.plot(t, measured_voltage, alpha=0.5, label='原始') plt.plot(t, filtered_voltage, 'r', label='陷波滤波后') plt.legend(); plt.title('时域信号对比') plt.subplot(2,1,2) plt.semilogy(f, Pxx, label='原始频谱') plt.semilogy(f_filt, Pxx_filt, 'r', label='滤波后频谱') plt.axvline(60, color='k', linestyle='--', label='60Hz干扰') plt.legend(); plt.xlabel('频率(Hz)'); plt.ylabel('功率谱密度') plt.tight_layout()

4.3 重采样与移动平均组合方案

当采样率与干扰频率存在整数倍关系时，重采样技术能显著提升滤波效果：

# 计算最优重采样率 target_fs = 17 * 60 # 1020Hz resampled_voltage = signal.resample(measured_voltage, int(len(measured_voltage)*target_fs/fs)) # 应用移动平均滤波 window_size = 17 smoothed_voltage = np.convolve( resampled_voltage, np.ones(window_size)/window_size, mode='same') # 可视化对比 t_resampled = np.arange(len(resampled_voltage)) / target_fs plt.figure(figsize=(12,4)) plt.plot(t, measured_voltage, alpha=0.3, label='原始信号') plt.plot(t_resampled, smoothed_voltage, 'r', label='重采样+移动平均') plt.legend(); plt.title('重采样技术对60Hz噪声���抑制效果')

5. 工程实践中的性能优化

5.1 实时处理与延迟补偿

在实际IoT应用中，我们需要考虑实时处理的计算效率：

from collections import deque class RealTimeFilter: def __init__(self, window_size=24): self.window = deque(maxlen=window_size) self.sum = 0.0 def update(self, new_value): if len(self.window) == self.window.maxlen: self.sum -= self.window[0] self.window.append(new_value) self.sum += new_value return self.sum / len(self.window) # 模拟实时处理 rt_filter = RealTimeFilter(window_size=24) online_filtered = [rt_filter.update(v) for v in tempC] plt.figure(figsize=(12,6)) plt.plot(temp_data.index, tempC, alpha=0.3, label='原始数据') plt.plot(temp_data.index, online_filtered, 'r', label='实时移动平均') plt.legend(); plt.title('实时滤波实现方案')

5.2 多传感器数据融合

当有多个温度传感器时，数据融合能提升结果可靠性：

# 模拟三个传感器的数据 sensors = { 'sensor1': tempC + np.random.normal(0, 0.5, len(tempC)), 'sensor2': tempC + np.random.normal(0, 0.8, len(tempC)), 'sensor3': tempC + np.random.normal(0, 1.0, len(tempC)) } # 卡尔曼滤波实现多传感器融合 from pykalman import KalmanFilter kf = KalmanFilter(transition_matrices=[1], observation_matrices=[1], initial_state_mean=tempC[0], observation_covariance=1.0, transition_covariance=0.1) state_means, _ = kf.filter(np.array(list(sensors.values())).T) plt.figure(figsize=(12,6)) for name, values in sensors.items(): plt.plot(temp_data.index, values, alpha=0.2, label=name) plt.plot(temp_data.index, state_means, 'k', linewidth=2, label='融合结果') plt.legend(); plt.title('多传感器数据融合效果')

查看全文

http://www.jsqmd.com/news/913773/