力传感器信号滤波处理方法
力传感器信号滤波处理方法。力传感器信号通常包含噪声、振动干扰和漂移等问题，需要针对性的滤波技术。
1. 力信号特性分析与滤波方案选择
function force_filter_demo()
% 力传感器信号滤波演示
% 生成模拟力信号
[t, force_clean, force_noisy] = generate_force_signal();
% 分析信号特性
analyze_force_signal(t, force_noisy);
% 应用多种滤波方法
filters = apply_all_filters(t, force_noisy);
% 比较滤波效果
compare_filter_performance(t, force_clean, force_noisy, filters);
% 实时滤波模拟
real_time_filter_simulation();
end
function [t, force_clean, force_noisy] = generate_force_signal()
% 生成模拟力信号：包含趋势、周期成分和冲击
Fs = 1000; % 采样频率 1kHz
t = 0:1/Fs:10; % 10秒数据
N = length(t);
% 1. 基础力信号（趋势+周期性）
trend = 100 + 5 * sin(2*pi*0.1*t); % 缓慢变化的趋势
periodic = 20 * sin(2*pi*5*t) + 10 * sin(2*pi*12*t); % 工作周期
% 2. 冲击信号（机械冲击）
impact_times = [2.5, 5.8, 8.2];
impacts = zeros(size(t));
for i = 1:length(impact_times)
idx = find(t >= impact_times(i), 1);
if ~isempty(idx)
% 高斯脉冲模拟冲击
pulse = 50 * exp(-((t - impact_times(i))/0.01).^2);
impacts = impacts + pulse;
end
end
% 3. 高频噪声（测量噪声）
hf_noise = 3 * randn(size(t));
% 4. 低频漂移（温度漂移）
drift = 2 * cumsum(randn(size(t))) / Fs;
% 合成信号
force_clean = trend + periodic + impacts;
force_noisy = force_clean + hf_noise + drift;
fprintf('生成力信号完成\n');
fprintf('信号长度: %d\n', N);
fprintf('采样频率: %.1f Hz\n', Fs);
end
function analyze_force_signal(t, force_signal)
% 分析力信号特性
Fs = 1/(t(2)-t(1));
N = length(force_signal);
figure('Position', [100, 100, 1200, 800]);
% 时域分析
subplot(2, 3, 1);
plot(t, force_signal, 'b');
title('原始力信号 - 时域');
xlabel('时间 (s)');
ylabel('力 (N)');
grid on;
% 频域分析
subplot(2, 3, 2);
[f, P] = compute_spectrum(force_signal, Fs);
semilogx(f, 10*log10(P), 'r');
title('功率谱密度');
xlabel('频率 (Hz)');
ylabel('功率 (dB)');
grid on;
% 统计特性
subplot(2, 3, 3);
histogram(force_signal, 50, 'FaceColor', 'green', 'FaceAlpha', 0.7);
title('信号分布');
xlabel('力值 (N)');
ylabel('频次');
grid on;
% 自相关分析
subplot(2, 3, 4);
[acf, lags] = xcorr(force_signal - mean(force_signal), 200, 'coeff');
plot(lags/Fs, acf, 'm');
title('自相关函数');
xlabel('时延 (s)');
ylabel('自相关系数');
grid on;
% 时频分析 - 短时傅里叶变换
subplot(2, 3, 5);
window = hamming(256);
noverlap = 200;
nfft = 1024;
spectrogram(force_signal, window, noverlap, nfft, Fs, 'yaxis');
title('时频分析 - 谱图');
% 信号特征统计
subplot(2, 3, 6);
stats = compute_signal_statistics(force_signal);
text(0.1, 0.9, '信号统计特征:', 'FontSize', 12, 'FontWeight', 'bold');
text(0.1, 0.7, sprintf('均值: %.2f N', stats.mean));
text(0.1, 0.6, sprintf('标准差: %.2f N', stats.std));
text(0.1, 0.5, sprintf('峰值: %.2f N', stats.peak));
text(0.1, 0.4, sprintf('峰峰值: %.2f N', stats.peak_to_peak));
text(0.1, 0.3, sprintf('信噪比: %.2f dB', stats.snr));
xlim([0, 1]);
ylim([0, 1]);
set(gca, 'XTick', [], 'YTick', []);
box on;
fprintf('\n信号分析完成\n');
end
2. 多种滤波方法实现
function filters = apply_all_filters(t, force_noisy)
% 应用多种滤波方法
Fs = 1/(t(2)-t(1));
filters = struct();
% 1. 移动平均滤波
filters.moving_avg = moving_average_filter(force_noisy, 20);
% 2. 中值滤波（对脉冲噪声有效）
filters.median = medfilt1(force_noisy, 15);
% 3. 巴特沃斯低通滤波
filters.butterworth_lowpass = butterworth_lowpass(force_noisy, Fs, 30);
% 4. 卡尔曼滤波
filters.kalman = kalman_filter(force_noisy);
% 5. 小波去噪
filters.wavelet = wavelet_denoise(force_noisy);
% 6. Savitzky-Golay滤波
filters.sgolay = sgolay_filter(force_noisy, 5, 21);
% 7. 自适应滤波
filters.adaptive = adaptive_lms_filter(force_noisy, Fs);
% 8. 零相位滤波（双向滤波）
filters.zero_phase = zero_phase_filter(force_noisy, Fs, 30);
fprintf('所有滤波方法应用完成\n');
end
function filtered = moving_average_filter(signal, window_size)
% 移动平均滤波
b = (1/window_size) * ones(1, window_size);
a = 1;
filtered = filter(b, a, signal);
% 处理边界效应
filtered(1:window_size-1) = signal(1:window_size-1);
end
function filtered = butterworth_lowpass(signal, Fs, cutoff_freq)
% 巴特沃斯低通滤波
order = 4; % 滤波器阶数
nyquist = Fs / 2;
Wn = cutoff_freq / nyquist;
[b, a] = butter(order, Wn, 'low');
filtered = filtfilt(b, a, signal); % 零相位滤波
end
function filtered = kalman_filter(signal)
% 简化卡尔曼滤波
% 状态变量: [位置; 速度]
% 观测变量: 力值
n = length(signal);
filtered = zeros(size(signal));
% 初始状态
x = [signal(1); 0]; % [位置; 速度]
P = [1, 0; 0, 1]; % 状态协方差
% 系统模型
A = [1, 1; 0, 1]; % 状态转移矩阵
H = [1, 0]; % 观测矩阵
Q = [0.1, 0; 0, 0.1]; % 过程噪声
R = 1; % 观测噪声
for k = 1:n
% 预测步骤
x = A * x;
P = A * P * A' + Q;
% 更新步骤
K = P * H' / (H * P * H' + R);
x = x + K * (signal(k) - H * x);
P = (eye(2) - K * H) * P;
filtered(k) = x(1);
end
end
function denoised = wavelet_denoise(signal)
% 小波去噪
% 使用小波变换分离噪声和信号
% 选择小波基
wavelet = 'db4';
level = 5;
% 小波分解
[C, L] = wavedec(signal, level, wavelet);
% 使用软阈值去噪
sigma = median(abs(C)) / 0.6745;
threshold = sigma * sqrt(2 * log(length(signal)));
% 应用阈值
C_denoised = wthresh(C, 's', threshold);
% 小波重构
denoised = waverec(C_denoised, L, wavelet);
% 确保长度一致
if length(denoised) > length(signal)
denoised = denoised(1:length(signal));
elseif length(denoised) < length(signal)
denoised(end+1:length(signal)) = denoised(end);
end
end
function filtered = sgolay_filter(signal, order, frame_size)
% Savitzky-Golay滤波（多项式平滑）
% 对峰值保持较好
if mod(frame_size, 2) == 0
frame_size = frame_size + 1; % 确保为奇数
end
filtered = sgolayfilt(signal, order, frame_size);
end
function filtered = adaptive_lms_filter(signal, Fs)
% 自适应LMS滤波
% 适用于时变噪声特性
n = length(signal);
mu = 0.01; % 步长
order = 10; % 滤波器阶数
% 生成参考噪声（这里用高频分量模拟）
d = signal;
x = sin(2*pi*60*t(1:n)) + 0.5*randn(1,n); % 参考噪声
% LMS滤波
w = zeros(order, 1); % 滤波器权重
filtered = zeros(size(signal));
for k = order:n
x_vec = x(k:-1:k-order+1);
y = w' * x_vec';
e = d(k) - y;
w = w + mu * e * x_vec';
filtered(k) = y;
end
% 处理前缘
filtered(1:order-1) = signal(1:order-1);
end
function filtered = zero_phase_filter(signal, Fs, cutoff_freq)
% 零相位滤波（双向滤波，无相位失真）
order = 4;
nyquist = Fs / 2;
Wn = cutoff_freq / nyquist;
[b, a] = butter(order, Wn, 'low');
filtered = filtfilt(b, a, signal);
end
3. 滤波性能评估与比较
function compare_filter_performance(t, force_clean, force_noisy, filters)
% 比较不同滤波方法的性能
filter_names = fieldnames(filters);
n_filters = length(filter_names);
% 计算各滤波器的性能指标
performance = struct();
for i = 1:n_filters
filter_name = filter_names{i};
filtered_signal = filters.(filter_name);
% 计算性能指标
perf = evaluate_filter_performance(force_clean, force_noisy, filtered_signal);
performance.(filter_name) = perf;
end
% 绘制比较结果
plot_filter_comparison(t, force_clean, force_noisy, filters, performance);
% 显示性能表格
display_performance_table(filter_names, performance);
end
function perf = evaluate_filter_performance(clean, noisy, filtered)
% 评估滤波性能
% 均方根误差
perf.rmse = sqrt(mean((filtered - clean).^2));
% 信噪比改善
snr_original = 10 * log10(var(clean) / var(noisy - clean));
snr_filtered = 10 * log10(var(clean) / var(filtered - clean));
perf.snr_improvement = snr_filtered - snr_original;
% 相关系数
perf.correlation = corr(clean(:), filtered(:));
% 峰值保持误差
[clean_peaks, clean_locs] = findpeaks(clean, 'MinPeakHeight', mean(clean));
[filtered_peaks, filtered_locs] = findpeaks(filtered, 'MinPeakHeight', mean(filtered));
if length(clean_peaks) > 0 && length(filtered_peaks) > 0
% 找到对应的峰值
peak_error = 0;
matched_peaks = 0;
for i = 1:length(clean_locs)
[~, idx] = min(abs(filtered_locs - clean_locs(i)));
if abs(filtered_locs(idx) - clean_locs(i)) < 10 % 在10个采样点内
peak_error = peak_error + abs(filtered_peaks(idx) - clean_peaks(i));
matched_peaks = matched_peaks + 1;
end
end
if matched_peaks > 0
perf.peak_error = peak_error / matched_peaks;
else
perf.peak_error = inf;
end
else
perf.peak_error = inf;
end
% 计算执行时间（这里用信号长度代替实际时间）
perf.computational_cost = length(filtered);
end
function plot_filter_comparison(t, force_clean, force_noisy, filters, performance)
% 绘制滤波结果比较
filter_names = fieldnames(filters);
n_filters = length(filter_names);
figure('Position', [100, 100, 1400, 1000]);
% 1. 原始信号与滤波信号对比
subplot(3, 3, 1);
plot(t, force_noisy, 'k', 'LineWidth', 0.5, 'DisplayName', '含噪信号');
hold on;
plot(t, force_clean, 'g', 'LineWidth', 2, 'DisplayName', '干净信号');
colors = lines(n_filters);
for i = 1:n_filters
filter_name = filter_names{i};
plot(t, filters.(filter_name), 'Color', colors(i,:), 'LineWidth', 1.5, ...
'DisplayName', sprintf('%s (RMSE=%.3f)', filter_name, performance.(filter_name).rmse));
end
title('滤波效果对比');
xlabel('时间 (s)');
ylabel('力 (N)');
legend('show', 'Location', 'best');
grid on;
% 2. 误差比较
subplot(3, 3, 2);
errors = zeros(n_filters, 1);
for i = 1:n_filters
filter_name = filter_names{i};
errors(i) = performance.(filter_name).rmse;
end
bar(errors);
set(gca, 'XTickLabel', filter_names);
ylabel('RMSE');
title('均方根误差比较');
grid on;
rotateXLabels(gca, 45);
% 3. 信噪比改善
subplot(3, 3, 3);
snr_imp = zeros(n_filters, 1);
for i = 1:n_filters
filter_name = filter_names{i};
snr_imp(i) = performance.(filter_name).snr_improvement;
end
bar(snr_imp);
set(gca, 'XTickLabel', filter_names);
ylabel('SNR改善 (dB)');
title('信噪比改善');
grid on;
rotateXLabels(gca, 45);
% 4. 相关系数
subplot(3, 3, 4);
corrs = zeros(n_filters, 1);
for i = 1:n_filters
filter_name = filter_names{i};
corrs(i) = performance.(filter_name).correlation;
end
bar(corrs);
set(gca, 'XTickLabel', filter_names);
ylabel('相关系数');
title('与真实信号相关性');
grid on;
rotateXLabels(gca, 45);
% 5. 峰值保持误差
subplot(3, 3, 5);
peak_errs = zeros(n_filters, 1);
for i = 1:n_filters
filter_name = filter_names{i};
if isfinite(performance.(filter_name).peak_error)
peak_errs(i) = performance.(filter_name).peak_error;
else
peak_errs(i) = 0;
end
end
bar(peak_errs);
set(gca, 'XTickLabel', filter_names);
ylabel('平均峰值误差 (N)');
title('峰值保持能力');
grid on;
rotateXLabels(gca, 45);
% 6. 频域比较
subplot(3, 3, 6);
Fs = 1/(t(2)-t(1));
[f_clean, P_clean] = compute_spectrum(force_clean, Fs);
plot(f_clean, 10*log10(P_clean), 'g', 'LineWidth', 2, 'DisplayName', '干净信号');
hold on;
for i = 1:n_filters
filter_name = filter_names{i};
[f_filt, P_filt] = compute_spectrum(filters.(filter_name), Fs);
plot(f_filt, 10*log10(P_filt), 'Color', colors(i,:), 'LineWidth', 1, ...
'DisplayName', filter_name);
end
xlabel('频率 (Hz)');
ylabel('功率 (dB)');
title('频域特性比较');
legend('show', 'Location', 'best');
grid on;
xlim([0, 100]);
% 7. 残差分析
subplot(3, 3, 7);
residuals = zeros(n_filters, length(t));
for i = 1:n_filters
filter_name = filter_names{i};
residuals(i, :) = filters.(filter_name) - force_clean;
end
boxplot(residuals', 'Labels', filter_names);
ylabel('残差 (N)');
title('滤波残差分布');
grid on;
rotateXLabels(gca, 45);
% 8. 计算复杂度比较
subplot(3, 3, 8);
costs = zeros(n_filters, 1);
for i = 1:n_filters
filter_name = filter_names{i};
costs(i) = performance.(filter_name).computational_cost;
end
bar(costs);
set(gca, 'XTickLabel', filter_names);
ylabel('相对计算成本');
title('计算复杂度');
grid on;
rotateXLabels(gca, 45);
% 9. 综合评分
subplot(3, 3, 9);
scores = zeros(n_filters, 1);
for i = 1:n_filters
filter_name = filter_names{i};
perf = performance.(filter_name);
% 综合评分公式（可根据需求调整权重）
score = (1/perf.rmse) * 0.4 + ... % RMSE权重40%
perf.snr_improvement * 0.2 + ... % SNR改善权重20%
perf.correlation * 0.2 + ... % 相关性权重20%
(1/(perf.peak_error+0.1)) * 0.1 + ... % 峰值保持权重10%
(1/(perf.computational_cost/1000)) * 0.1; % 计算成本权重10%
scores(i) = score;
end
% 归一化评分
scores = scores / max(scores) * 100;
bar(scores);
set(gca, 'XTickLabel', filter_names);
ylabel('综合评分 (%)');
title('滤波器综合性能评分');
grid on;
rotateXLabels(gca, 45);
% 添加数值标签
for i = 1:length(scores)
text(i, scores(i), sprintf('%.1f', scores(i)), ...
'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');
end
end
4. 实用工具函数
function [f, P] = compute_spectrum(signal, Fs)
% 计算信号功率谱
N = length(signal);
f = (0:N/2) * Fs / N;
Y = fft(signal);
P = abs(Y(1:N/2+1)).^2 / N;
end
function stats = compute_signal_statistics(signal)
% 计算信号统计特征
stats.mean = mean(signal);
stats.std = std(signal);
stats.peak = max(signal);
stats.peak_to_peak = max(signal) - min(signal);
% 估计信噪比（假设信号主要成分在低频）
[f, P] = compute_spectrum(signal, 1000);
signal_power = sum(P(f < 50)); % 50Hz以下为信号
noise_power = sum(P(f >= 50)); % 50Hz以上为噪声
stats.snr = 10 * log10(signal_power / noise_power);
end
function display_performance_table(filter_names, performance)
% 显示性能比较表格
fprintf('\n=== 滤波器性能比较 ===\n');
fprintf('%-20s %-8s %-12s %-12s %-12s\n', ...
'滤波器', 'RMSE', 'SNR改善(dB)', '相关系数', '峰值误差');
fprintf('%s\n', repmat('-', 1, 70));
for i = 1:length(filter_names)
filter_name = filter_names{i};
perf = performance.(filter_name);
if isfinite(perf.peak_error)
peak_err_str = sprintf('%.3f', perf.peak_error);
else
peak_err_str = 'NaN';
end
fprintf('%-20s %-8.3f %-12.3f %-12.3f %-12s\n', ...
filter_name, perf.rmse, perf.snr_improvement, ...
perf.correlation, peak_err_str);
end
end
function real_time_filter_simulation()
% 实时滤波模拟
fprintf('\n开始实时滤波模拟...\n');
Fs = 1000; % 采样频率
simulation_time = 5; % 秒
t_rt = 0:1/Fs:simulation_time;
% 生成实时信号
signal_rt = 100 + 10 * sin(2*pi*2*t_rt) + 5 * sin(2*pi*8*t_rt) + ...
3 * randn(size(t_rt)) + 0.5 * cumsum(randn(size(t_rt))) / Fs;
% 实时滤波参数
buffer_size = 100; % 缓冲区大小
filtered_rt = zeros(size(signal_rt));
figure('Position', [200, 200, 1000, 600]);
for i = buffer_size:length(signal_rt)
% 获取当前缓冲区
if i <= buffer_size
buffer = signal_rt(1:i);
else
buffer = signal_rt(i-buffer_size+1:i);
end
% 应用实时滤波（移动平均）
filtered_rt(i) = mean(buffer);
% 实时显示（每100个点更新一次）
if mod(i, 100) == 0
clf;
subplot(2,1,1);
plot(t_rt(1:i), signal_rt(1:i), 'b', 'LineWidth', 1, 'DisplayName', '原始信号');
hold on;
plot(t_rt(1:i), filtered_rt(1:i), 'r', 'LineWidth', 2, 'DisplayName', '滤波后');
title(sprintf('实时滤波 (时间: %.1f/%.1f s)', t_rt(i), simulation_time));
xlabel('时间 (s)');
ylabel('力 (N)');
legend('show');
grid on;
subplot(2,1,2);
recent_idx = max(1, i-500):i;
plot(t_rt(recent_idx), signal_rt(recent_idx), 'b', 'LineWidth', 1);
hold on;
plot(t_rt(recent_idx), filtered_rt(recent_idx), 'r', 'LineWidth', 2);
title('最近0.5秒数据');
xlabel('时间 (s)');
ylabel('力 (N)');
grid on;
drawnow;
end
end
fprintf('实时滤波模拟完成\n');
end
% 辅助函数：旋转X轴标签
function rotateXLabels(ax, angle)
if ~exist('angle', 'var')
angle = 45;
end
set(ax, 'XTickLabelRotation', angle);
end
5. 实际应用建议
function practical_filtering_advice()
% 提供实际滤波建议
fprintf('\n=== 力传感器滤波实用建议 ===\n');
advice = {
'1. 移动平均滤波: 适合平滑高频噪声，计算简单，实时性好'
'2. 中值滤波: 对脉冲噪声有效，能保持边缘'
'3. 巴特沃斯滤波: 频域控制精确，相位线性'
'4. 卡尔曼滤波: 适合时变系统，需要系统模型'
'5. 小波去噪: 多尺度分析，适合非平稳信号'
'6. Savitzky-Golay: 峰值保持好，适合光谱数据'
'7. 自适应滤波: 适合噪声特性变化的场景'
'8. 零相位滤波: 无相位失真，适合时序分析'
''
'选择原则:'
'- 实时应用: 移动平均、中值滤波'
'- 精确控制: 巴特沃斯、零相位滤波'
'- 冲击信号: 中值滤波、小波去噪'
'- 漂移消除: 高通滤波、自适应滤波'
'- 综合性能: 卡尔曼滤波、小波去噪'
};
for i = 1:length(advice)
fprintf('%s\n', advice{i});
end
end
% 运行完整演示
if __name__ == "__main__"
force_filter_demo();
practical_filtering_advice();
end
参考代码 对力传感器测量到的力信号进行滤波 www.youwenfan.com/contentcni/65779.html
总结
- 全面的滤波方法：包含8种常用滤波算法
- 性能评估体系：多维度评估滤波效果
- 实时处理能力：支持实时滤波模拟
- 可视化分析：完整的信号分析和结果展示
- 实用建议：根据应用场景推荐合适方法
浙公网安备 33010602011771号