一、算法原理概述

1. 小波變換（Wavelet Transform）

   • 通過多尺度分解將信號分為高頻（細節）和低頻（近似）部分，高頻通常包含噪聲，低頻保留主體信息。
   • 使用Haar小波（計算高效）進行快速分解與重構：

     // Haar小波分解公式
     approx = (x[2i] + x[2i+1]) / 2.0;
     detail = (x[2i] - x[2i+1]) / 2.0;
     

2. 卡爾曼濾波（Kalman Filter）

   • 預測-更新兩階段遞歸估計[citation:5][citation:6]：
     • 預測：基于狀態方程預估當前狀態和誤差協方差

       Xpred=A?Xprev+B?uPpred=A?Pprev?AT+QXpred?=A?Xprev?+B?uPpred?=A?Pprev??AT+Q

     • 更新：結合觀測值修正預測值，計算卡爾曼增益K

       K=Ppred?HT/(H?Ppred?HT+R)Xnew=Xpred+K?(z?H?Xpred)Pnew=(I?K?H)?PpredK=Ppred??HT/(H?Ppred??HT+R)Xnew?=Xpred?+K?(z?H?Xpred?)Pnew?=(I?K?H)?Ppred?

3. 小波卡爾曼融合

   • 先對原始信號小波分解，對不同頻段獨立應用卡爾曼濾波，最后重構信號
   • 高頻部分使用更高測量噪聲協方差R（抑制噪聲），低頻部分使用更低R（保留主體）。

二、Rust實現代碼

cargo.toml

[package]
name = "wavelet-kalman-filter"
version = "0.1.0"
edition = "2021"[dependencies]
plotters = "0.3.5"
rand = "0.8.5"
log = "0.4.17"
simple_logger = "5.0.0"

main.rs

use log::{info};
use log::LevelFilter;
use simple_logger::SimpleLogger;// 初始化日志記錄，將日志保存到文件
fn init_logger() {SimpleLogger::new().with_level(LevelFilter::Info).with_utc_timestamps().init().expect("Failed to initialize logger");info!("日志系統初始化完成");
}// 主結構體
pub struct WaveletKalmanFilter {pub a: f64,      // 狀態轉移系數（如勻速運動時A=1）pub q: f64,      // 過程噪聲協方差pub r_low: f64,  // 低頻觀測噪聲協方差pub r_high: f64, // 高頻觀測噪聲協方差pub state: f64,  // 當前狀態估計pub cov: f64,    // 當前誤差協方差
}// 一維Haar小波分解（返回低頻近似 + 高頻細節）
pub fn haar_decompose(signal: &[f64]) -> (Vec<f64>, Vec<f64>) {info!("開始Haar小波分解，輸入信號長度: {}", signal.len());// 對信號進行對稱延拓，以緩解小波分析的邊界問題fn symmetric_extend(signal: &[f64]) -> Vec<f64> {if signal.len() <= 1 {info!("信號長度小于等于1，直接返回原信號");return signal.to_vec();}let mut extended = Vec::new();// 左側對稱部分for i in (1..=signal.len() - 1).rev() {extended.push(signal[i]);}// 原始信號部分extended.extend_from_slice(signal);// 右側對稱部分for i in 1..signal.len() {extended.push(signal[signal.len() - 1 - i]);}info!("對稱延拓后信號長度: {} 信號內容: {:?}", extended.len(), extended);extended}// 對延拓后的信號進行Haar小波分解pub fn haar_decompose(signal: &[f64]) -> (Vec<f64>, Vec<f64>) {let extended = symmetric_extend(signal);let mut approx = Vec::new();let mut detail = Vec::new();for i in (0..extended.len()).step_by(2) {let sum = extended[i] + extended.get(i + 1).unwrap_or(&0.0);let diff = extended[i] - extended.get(i + 1).unwrap_or(&0.0);approx.push(sum / 2.0);detail.push(diff / 2.0);}// 截取與原始信號分解后相同長度的結果let half_len = (signal.len() + 1) / 2;let approx = approx[extended.len() / 4..extended.len() / 4 + half_len].to_vec();let detail = detail[extended.len() / 4..extended.len() / 4 + half_len].to_vec();info!("延拓后Haar分解完成，近似分量長度: {}, 細節分量長度: {}",approx.len(),detail.len());(approx, detail)}let mut approx = Vec::new();let mut detail = Vec::new();for i in (0..signal.len()).step_by(2) {let sum = signal[i] + signal.get(i + 1).unwrap_or(&0.0);let diff = signal[i] - signal.get(i + 1).unwrap_or(&0.0);approx.push(sum / 2.0);detail.push(diff / 2.0);}info!("直接Haar分解完成，近似分量長度: {}, 細節分量長度: {}",approx.len(),detail.len());(approx, detail)
}// 一維Haar小波重構
pub fn haar_recompose(approx: &[f64], detail: &[f64]) -> Vec<f64> {info!("開始Haar小波重構，近似分量長度: {}, 細節分量長度: {}",approx.len(),detail.len());let mut output = Vec::with_capacity(approx.len() * 2);for i in 0..approx.len() {output.push(approx[i] + detail[i]);output.push(approx[i] - detail[i]);}info!("Haar小波重構完成，輸出信號長度: {}", output.len());output
}impl WaveletKalmanFilter {// 初始化濾波器pub fn new(a: f64, q: f64, r_low: f64, r_high: f64) -> Self {info!("創建新的WaveletKalmanFilter，a: {}, q: {}, r_low: {}, r_high: {}",a, q, r_low, r_high);WaveletKalmanFilter {a,q,r_low,r_high,state: 0.0,cov: 1.0, // 初始協方差設為較大值（表示不確定性高）}}// 單步預測pub fn predict(&mut self) {info!("執行單步預測前，狀態: {}, 協方差: {}", self.state, self.cov);self.state = self.a * self.state; // X_pred = A * X_prevself.cov = self.a * self.cov * self.a + self.q; // P_pred = A*P*A^T + Qinfo!("執行單步預測后，狀態: {}, 協方差: {}", self.state, self.cov);}// 單步更新（根據頻段選擇R）pub fn update(&mut self, measurement: f64, is_low_freq: bool) {// info!(//     "執行單步更新前，測量值: {}, 是否低頻: {}",//     measurement, is_low_freq// );let r = if is_low_freq { self.r_low } else { self.r_high };let k = self.cov / (self.cov + r); // 卡爾曼增益Kself.state += k * (measurement - self.state); // X_new = X_pred + K*(z - H*X_pred)self.cov = (1.0 - k) * self.cov; // P_new = (I - K*H)*P_predinfo!("執行單步更新后，狀態: {}, 協方差: {}", self.state, self.cov);}
}pub fn wavelet_kalman_filter(signal: &[f64],levels: usize, // 小波分解層數a: f64,        // 狀態轉移系數q: f64,        // 過程噪聲r_low: f64,    // 低頻觀測噪聲r_high: f64,   // 高頻觀測噪聲
) -> Vec<f64> {info!("開始小波卡爾曼濾波，輸入信號長度: {}, 分解層數: {}",signal.len(),levels);// 1. 遞歸小波分解（多尺度）let mut approx = signal.to_vec();let mut details = Vec::new();for level in 0..levels {info!("第 {} 層小波分解", level + 1);let (a, d) = haar_decompose(&approx);details.push(d);approx = a;}// 2. 對各頻段獨立應用卡爾曼濾波let mut kalman = WaveletKalmanFilter::new(a, q, r_low, r_high);kalman.state = approx[0]; // 初始化狀態// 低頻濾波（R較小，信任觀測值）info!("開始低頻濾波，近似分量長度: {}", approx.len());for i in 0..approx.len() {info!("低頻濾波第 {} 步", i + 1);kalman.predict();kalman.update(approx[i], true); // 標記為低頻approx[i] = kalman.state;}// 高頻濾波（R較大，抑制噪聲）info!("開始高頻濾波，細節分量數量: {}", details.len());for (level, d) in details.iter_mut().enumerate() {info!("第 {} 層細節分量濾波，長度: {}", level + 1, d.len());kalman = WaveletKalmanFilter::new(a, q, r_low, r_high); // 重置濾波器kalman.state = d[0];for i in 0..d.len() {info!("第 {} 層細節分量第 {} 步濾波", level + 1, i + 1);kalman.predict();kalman.update(d[i], false); // 標記為高頻d[i] = kalman.state;}}// 3. 小波重構info!("開始小波重構，分解層數: {}", levels);for _ in 0..levels {approx = haar_recompose(&approx, &details.pop().unwrap());}info!("小波卡爾曼濾波完成，輸出信號長度: {}", approx.len());approx
}// 繪圖函數實現
fn draw_signals(clean: &[f64],noisy: &[f64],filtered: &[f64],
) -> Result<(), Box<dyn std::error::Error>> {info!("開始繪制信號圖，原始信號長度: {}, 含噪信號長度: {}, 濾波信號長度: {}",clean.len(),noisy.len(),filtered.len());use plotters::prelude::*;let root = BitMapBackend::new("signals.png", (1024, 768)).into_drawing_area();root.fill(&WHITE)?;let mut chart = ChartBuilder::on(&root).caption("信號處理結果", ("sans-serif", 50).into_font()).margin(10).x_label_area_size(30).y_label_area_size(30).build_cartesian_2d(0f64..clean.len() as f64, -2f64..2f64)?;chart.configure_mesh().draw()?;chart.draw_series(LineSeries::new(clean.iter().enumerate().map(|(i, &y)| (i as f64, y)),&RED,))?.label("原始信號").legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &RED));chart.draw_series(LineSeries::new(noisy.iter().enumerate().map(|(i, &y)| (i as f64, y)),&BLUE,))?.label("含噪信號").legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &BLUE));chart.draw_series(LineSeries::new(filtered.iter().enumerate().map(|(i, &y)| (i as f64, y)),&GREEN,))?.label("濾波信號").legend(|(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], &GREEN));chart.configure_series_labels().background_style(&WHITE.mix(0.8)).border_style(&BLACK).draw()?;info!("信號圖繪制完成，保存至 signals.png");Ok(())
}// 測試函數：添加高斯噪聲的信號濾波
fn test_noisy_sine() {info!("開始測試含噪正弦信號濾波");let clean_signal: Vec<_> = (0..1000).map(|i| (i as f64 * 0.1).sin()).collect();let noisy_signal: Vec<_> = clean_signal.iter().map(|x| x + rand::random::<f64>()).collect();let filtered = wavelet_kalman_filter(&noisy_signal, 3, 1.0, 0.01, 0.1, 2.0);// 繪圖比較（使用plotters庫）draw_signals(&clean_signal, &noisy_signal, &filtered);info!("含噪正弦信號濾波測試完成");
}// 主函數
fn main() {init_logger();info!("程序啟動");test_noisy_sine();info!("程序結束");
}

?

?