基于小波域隐马尔科夫模型的信号处理：从EM算法到MATLAB仿真实践

1. 项目缘起从“盲人摸象”到“庖丁解牛”的信号处理困境在信号处理的世界里我们常常扮演着“盲人摸象”的角色。面对一段复杂的信号比如一段夹杂着噪声的语音、一幅纹理丰富的医学影像或者一组非平稳的金融时间序列我们手头只有观测到的数据。这些数据背后信号的真实状态比如某个时刻是清音还是浊音、图像某个区域是边缘还是平滑、市场处于牛市还是熊市是隐藏的、不可直接观测的。更棘手的是这些隐藏的状态之间并非独立它们按照某种概率规律相互转移而每个隐藏状态又会以某种概率分布“发射”出我们观测到的数据。这就像观察一个人的行为观测去推测他的情绪隐藏状态而情绪会变化且同一种情绪可能表现出不同的行为。传统的单一模型比如直接用高斯模型去拟合整个信号或者用小波变换后简单阈值去噪往往力不从心。它们要么忽略了状态间的时序关联要么对信号局部特性的刻画过于粗糙。这就引出了我们这次要折腾的“组合拳”基于小波域隐马尔科夫模型Wavelet-Domain Hidden Markov Model, WD-HMM的参数估计。简单来说这个模型的思路堪称“庖丁解牛”。首先用小波变换这把“解剖刀”将信号从时域/空域变换到小波域。小波变换的好处是它能同时在时间和频率或空间和尺度上提供良好的局部化特性信号的特征如边缘、瞬态成分会被集中到少数几个系数上。然后我们对这些小波系数建立隐马尔科夫模型HMM。为什么是HMM因为小波系数之间存在跨尺度的统计依赖性——一个粗尺度低频的大系数很可能预示着其对应的细尺度高频子带上也会出现大系数这种依赖关系用马尔科夫链来描述非常合适。而“隐”则是因为我们假设每个小波系数属于某种“状态”例如“大系数”、“小系数”、“正边缘”、“负边缘”等这个状态是隐藏的但我们能观测到系数值。模型建好了核心问题来了模型的参数状态转移概率、状态初始分布、每个状态对应的高斯分布的均值和方差我们并不知道。这就需要期望最大化EM算法出场了。EM算法是解决这类含有隐变量模型参数估计问题的利器。它通过迭代执行两步E步期望步基于当前参数估计隐藏状态的后验概率M步最大化步利用E步得到的后验概率更新模型参数以最大化数据的期望似然。如此循环直至参数收敛。在MATLAB里仿真这套流程意义重大。它不仅是理论到实践的桥梁更能让我们直观地理解模型如何工作、参数如何迭代更新、以及最终模型对信号如去噪、分割、分类的提升效果。下面我就把自己在MATLAB中实现WD-HMM参数估计仿真时趟过的路、踩过的坑以及核心的思考逻辑毫无保留地分享出来。2. 模型深潜WD-HMM的数学骨架与EM算法的迭代舞步在动手写代码之前我们必须把模型的数学骨架搭清楚知道每一步计算究竟在干什么。否则代码就会变成一堆看不懂的魔法数字。2.1 小波域隐马尔科夫模型WD-HMM的形式化定义我们的观测数据是经过小波分解后得到的一组系数。假设我们对一维信号进行J层离散小波变换DWT会得到一个近似系数最粗尺度和J组细节系数。通常HMM建模主要针对细节系数因为近似系数包含了信号的概貌其统计特性不同。设第j尺度j1是最细尺度jJ是最粗尺度上有N_j个细节系数记为观测序列O {o_1, o_2, ..., o_T}这里T是所有尺度细节系数的总数按特定扫描顺序如从粗到细、同尺度内从左到右。实际上为了利用跨尺度依赖性我们常按小波树或小波四叉树对图像来组织系数和建立HMM。一个经典的HMM由以下五元组λ (S, V, A, B, π)定义在WD-HMM语境下S: 隐藏状态的集合。假设有K个状态S {s1, s2, ..., sK}。在信号处理中通常K2大/小系数或K3大正值/小值/大负值。V: 观测值的集合。在我们的连续观测HMM中观测是实数值的小波系数所以V是实数集R。这意味着B不再是离散概率矩阵而是概率密度函数。A: 状态转移概率矩阵。A [a_{ij}]其中a_{ij} P(q_{t1} s_j | q_t s_i)表示在时刻t处于状态s_i的条件下t1时刻转移到状态s_j的概率。在跨尺度HMM中“时刻t”对应着小波树中父节点的位置“时刻t1”可能对应其子节点位置转移概率刻画了状态从父节点到子节点的遗传特性。B: 观测概率密度函数。对于连续HMM通常假设给定状态s_k下观测值o_t服从高斯分布混合高斯更通用但更复杂即b_k(o_t) N(o_t; μ_k, σ_k^2)。其中μ_k和σ_k^2是状态s_k对应的高斯分布的均值和方差。π: 初始状态概率分布。π [π_i]其中π_i P(q_1 s_i)表示在序列起始时刻最粗尺度的根节点处于状态s_i的概率。我们的目标就是给定观测序列O估计出模型参数λ (A, B, π)其中B包含{μ_k, σ_k^2}。2.2 EM算法Baum-Welch算法的详细推演对于HMMEM算法的具体实现就是著名的Baum-Welch算法。它通过定义两个核心变量来进行迭代前向变量 α_t(i): 在给定模型λ下到时刻t为止的观测序列为(o_1, o_2, ..., o_t)且时刻t隐藏状态为s_i的概率。 α_t(i) P(o_1, o_2, ..., o_t, q_t s_i | λ) 初始化α_1(i) π_i * b_i(o_1) 递推α_{t1}(j) [Σ_{i1}^{K} α_t(i) * a_{ij}] * b_j(o_{t1}), for t1,2,...,T-1后向变量 β_t(i): 在给定模型λ和时刻t隐藏状态为s_i的条件下从t1到T的观测序列为(o_{t1}, ..., o_T)的概率。 β_t(i) P(o_{t1}, o_{t2}, ..., o_T | q_t s_i, λ) 初始化β_T(i) 1 (通常约定) 递推β_t(i) Σ_{j1}^{K} a_{ij} * b_j(o_{t1}) * β_{t1}(j), for tT-1, T-2, ..., 1基于α和β我们可以计算两个至关重要的期望统计量γ_t(i): 在给定模型λ和全部观测O的条件下时刻t隐藏状态为s_i的概率。 γ_t(i) P(q_t s_i | O, λ) [α_t(i) * β_t(i)] / Σ_{j1}^{K} α_t(j) * β_t(j) 这个量是E步的核心产出它是对隐藏状态的“软分配”。ξ_t(i, j): 在给定模型λ和全部观测O的条件下时刻t隐藏状态为s_i且时刻t1隐藏状态为s_j的概率。 ξ_t(i, j) P(q_t s_i, q_{t1} s_j | O, λ) [α_t(i) * a_{ij} * b_j(o_{t1}) * β_{t1}(j)] / P(O|λ) 其中P(O|λ) Σ_{i1}^{K} α_T(i)是整个观测序列的似然。有了γ_t(i)和ξ_t(i, j)M步的参数更新公式就非常直观了它们实际上是最大似然估计MLE在已知状态后验分布下的加权版本更新初始分布 π_iπ_i_new γ_1(i) 在序列开始时刻处于状态i的期望概率更新转移概率 a_{ij}a_{ij}new Σ{t1}^{T-1} ξ_t(i, j) / Σ_{t1}^{T-1} γ_t(i) 从状态i转移到j的期望次数除以离开状态i的总期望次数更新观测分布参数高斯均值 μ_k_newμ_k_new Σ_{t1}^{T} γ_t(k) * o_t / Σ_{t1}^{T} γ_t(k) 所有观测的加权平均权重是该观测属于状态k的概率方差 σ_k^2_newσ_k^2_new Σ_{t1}^{T} γ_t(k) * (o_t - μ_k_new)^2 / Σ_{t1}^{T} γ_t(k) 加权方差算法流程就是初始化参数λ^0 → (E步) 计算α, β, γ, ξ → (M步) 用上述公式更新参数得到λ^1 → 重复E步和M步直到对数似然log P(O|λ)的变化小于某个阈值或达到最大迭代次数。注意在计算α_t(i)时由于连乘很多小于1的概率其值会迅速下溢到0。必须使用缩放技巧Scaling。通常对每个时刻t的α_t(i)进行缩放令其和为1缩放因子记为c_t 1 / Σ_i α_t(i)。相应地β_t(i)也需要用相同的缩放因子进行补偿。最终序列的对数似然可以通过缩放因子计算log P(O|λ) -Σ_{t1}^{T} log(c_t)。忽略这个细节你的仿真会在几次迭代后就因数值下溢而崩溃。3. MATLAB仿真实战从信号生成到模型评估理论清晰后我们进入MATLAB实战环节。我将以一个一维非平稳信号例如一个阶跃信号加高斯白噪声的去噪为例展示完整的仿真流程。3.1 仿真环境搭建与测试信号生成首先我们生成一个包含突变的测试信号。这里不使用现成的噪声图像而是自己构造以便精确评估去噪效果。% 参数设置 clear; close all; clc; rng(42); % 固定随机种子确保结果可复现 N 1024; % 信号长度 t (0:N-1)/N; % 生成原始干净信号一个阶跃 一个脉冲 x_clean zeros(N,1); x_clean(300:500) 2; % 阶跃 x_clean(700) 3; % 脉冲 % 添加高斯白噪声 noise_std 0.5; % 噪声标准差 x_noisy x_clean noise_std * randn(N,1); % 可视化 figure; subplot(2,1,1); plot(t, x_clean); title(原始干净信号); grid on; ylim([-1 4]); subplot(2,1,2); plot(t, x_noisy); title([加噪信号 (噪声\sigma, num2str(noise_std), )]); grid on; ylim([-1 4]);3.2 小波变换与系数组织我们使用MATLAB的wavedec函数进行离散小波分解。选择合适的小波基和分解层数很重要。对于突变信号db4Daubechies 4或sym4Symlets 4是不错的选择它们在光滑性和紧支撑性之间有较好平衡。分解层数J通常选3-5层这里选4层。% 小波分解 wavelet_name db4; J 4; % 分解层数 [C, L] wavedec(x_noisy, J, wavelet_name); % C: 系数向量L: 各层系数长度记录 % 提取各层细节系数 detail_coeffs cell(1, J); for j 1:J detail_coeffs{j} detcoef(C, L, j); % 第j层细节系数 end approx_coeff appcoef(C, L, wavelet_name); % 近似系数 % 为了建立跨尺度HMM我们需要按小波树结构组织系数。 % 对于一维信号小波树是二叉树。我们将系数组织成一个长向量顺序是 % [第J层(最粗)细节系数, 第J层系数的第一个孩子(第J-1层对应位置), 第二个孩子, ...] % 这需要自己编写一个索引映射函数。这里为了简化演示我们采用一种近似 % 将所有细节系数按从粗尺度到细尺度的顺序拼接并假设一个简单的状态转移只与相邻尺度的“父-子”关系有关。 % 在实际高质量仿真中必须实现精确的树结构索引。 obs_sequence []; for j J:-1:1 % 从粗到细 obs_sequence [obs_sequence; detail_coeffs{j}(:)]; end T length(obs_sequence); % 观测序列长度实操心得wavedec返回的C和L结构是理解小波系数的关键。L是一个数组其元素依次是近似系数的长度第J层细节系数长度第J-1层细节系数长度...第1层细节系数长度。detcoef和appcoef函数依赖于这个L数组来正确提取系数。自己管理这些索引是后续构建树结构HMM的基础也是容易出错的地方。3.3 HMM参数初始化与Baum-Welch算法实现这是仿真的核心。我们假设有K2个隐藏状态状态1代表“小系数”对应噪声或平滑区域状态2代表“大系数”对应信号边缘或突变。K 2; % 隐藏状态数 max_iter 50; % EM最大迭代次数 tol 1e-6; % 收敛阈值 % 1. 初始化参数 % 初始分布假设开始时处于“小系数”状态的概率高 pi_init [0.8; 0.2]; % 转移矩阵对角元保持同一状态概率高非对角元切换状态概率低。 % 同时我们期望“大系数”状态更倾向于产生“大系数”子节点。 A_init [0.9, 0.1; 0.2, 0.8]; % 行和为1 % 观测高斯分布参数根据系数直方图粗略初始化 % 假设系数大致服从零均值高斯但“大系数”状态的方差更大可能包含信号。 coeff_std std(obs_sequence); coeff_mean mean(obs_sequence); % 状态1小均值接近0方差小接近噪声方差 mu_init [coeff_mean; coeff_mean 0.5*coeff_std]; sigma2_init [(0.5*coeff_std)^2; (1.5*coeff_std)^2]; % 将参数打包 lambda.A A_init; lambda.B.mu mu_init; lambda.B.sigma2 sigma2_init; % 存储方差 lambda.pi pi_init; % 2. 实现带缩放的Baum-Welch算法 log_likelihood_curve zeros(max_iter, 1); for iter 1:max_iter %% ---------- E步计算前向、后向变量及缩放因子 ---------- % 初始化 alpha zeros(T, K); beta zeros(T, K); scale_factor zeros(T, 1); % 缩放因子c_t % 计算观测概率密度高斯 B zeros(T, K); for k 1:K B(:, k) (1./sqrt(2*pi*lambda.B.sigma2(k))) .* exp(-0.5 * (obs_sequence - lambda.B.mu(k)).^2 ./ lambda.B.sigma2(k)); end % 避免零概率数值稳定 B(B eps) eps; % 前向算法 with scaling % t1 alpha(1, :) lambda.pi(:) .* B(1, :); scale_factor(1) 1 / sum(alpha(1, :)); alpha(1, :) alpha(1, :) * scale_factor(1); % t2:T for t 2:T for j 1:K alpha(t, j) sum(alpha(t-1, :) .* lambda.A(:, j)) * B(t, j); end scale_factor(t) 1 / sum(alpha(t, :)); alpha(t, :) alpha(t, :) * scale_factor(t); end % 计算当前迭代的对数似然 log_likelihood_curve(iter) -sum(log(scale_factor)); % 后向算法 with scaling (使用相同的scale_factor) beta(T, :) 1 * scale_factor(T); % 初始化并缩放 for t T-1:-1:1 for i 1:K beta(t, i) sum(lambda.A(i, :) .* B(t1, :) .* beta(t1, :)); end beta(t, :) beta(t, :) * scale_factor(t); end % 计算gamma_t(i) 和 xi_t(i,j) gamma zeros(T, K); for t 1:T gamma(t, :) alpha(t, :) .* beta(t, :); gamma(t, :) gamma(t, :) / sum(gamma(t, :)); % 归一化实际上由于缩放alpha.*beta已经归一化这里再加一次确保数值稳定 end xi zeros(T-1, K, K); for t 1:T-1 denom sum(sum( alpha(t, :) .* lambda.A .* (B(t1, :) .* beta(t1, :)) )); % P(O|λ)的缩放版本 for i 1:K for j 1:K xi(t, i, j) alpha(t, i) * lambda.A(i, j) * B(t1, j) * beta(t1, j) / denom; end end end %% ---------- M步重新估计参数 ---------- % 更新初始分布 lambda.pi gamma(1, :); % 更新转移矩阵 for i 1:K for j 1:K lambda.A(i, j) sum(squeeze(xi(:, i, j))) / sum(gamma(1:end-1, i)); end % 确保每行和为1数值计算可能导致微小偏差 lambda.A(i, :) lambda.A(i, :) / sum(lambda.A(i, :)); end % 更新高斯分布参数 for k 1:K gamma_k gamma(:, k); lambda.B.mu(k) sum(gamma_k .* obs_sequence) / sum(gamma_k); lambda.B.sigma2(k) sum(gamma_k .* (obs_sequence - lambda.B.mu(k)).^2) / sum(gamma_k); % 防止方差过小导致数值问题 lambda.B.sigma2(k) max(lambda.B.sigma2(k), eps); end %% ---------- 检查收敛 ---------- if iter 1 delta_loglik abs(log_likelihood_curve(iter) - log_likelihood_curve(iter-1)); if delta_loglik tol fprintf(EM算法在第 %d 次迭代收敛。\n, iter); log_likelihood_curve log_likelihood_curve(1:iter); % 截断 break; end end end if iter max_iter fprintf(达到最大迭代次数 %d。\n, max_iter); end % 绘制对数似然曲线 figure; plot(1:length(log_likelihood_curve), log_likelihood_curve, b-o, LineWidth, 1.5); xlabel(迭代次数); ylabel(对数似然); title(EM算法收敛曲线); grid on;3.4 状态解码与信号重构训练好模型后我们可以利用维特比Viterbi算法找到最可能的状态序列硬判决然后根据状态进行系数处理。例如对于去噪我们可以将属于“小系数”状态的系数置零或收缩而保留“大系数”状态的系数。% 维特比算法解码最可能状态序列 delta zeros(T, K); psi zeros(T, K, int32); % 记录回溯路径 % 初始化 delta(1, :) log(lambda.pi(:)) log(B(1, :)); psi(1, :) 0; % 递推 for t 2:T for j 1:K [max_val, max_idx] max(delta(t-1, :) log(lambda.A(:, j))); delta(t, j) max_val log(B(t, j)); psi(t, j) max_idx; end end % 终止与回溯 [~, best_last_state] max(delta(T, :)); best_path zeros(T, 1, int32); best_path(T) best_last_state; for t T-1:-1:1 best_path(t) psi(t1, best_path(t1)); end % 根据解码状态处理小波系数 % 假设状态1为“噪声/小系数”状态2为“信号/大系数” processed_coeffs_flat obs_sequence; processed_coeffs_flat(best_path 1) 0; % 将“小系数”状态置零硬阈值 % 或者使用软阈值processed_coeffs_flat(best_path 1) sign(obs_sequence(best_path1)) .* max(abs(obs_sequence(best_path1)) - threshold, 0); % 将处理后的平坦系数向量还原回小波分解结构C_recon % 这是仿真的另一个难点需要逆向进行3.2节中的系数组织操作。 % 这里我们简化处理直接修改原始系数向量C中对应细节系数的部分。 % 首先我们需要知道obs_sequence中每个元素对应原始C向量中的哪个位置。 % 这需要根据L数组和之前组织的顺序从粗到细来建立精确的索引映射。 % 由于篇幅和简化演示我们这里跳过精确映射假设一种简单情况。 % 在实际应用中必须实现这个逆向映射。 fprintf(注意此处系数逆向映射为简化演示。完整实现需精确的索引管理。\n); % 为了展示效果我们采用另一种更直接的思路利用训练好的模型计算每个系数的状态概率gamma % 然后对系数进行基于概率的软阈值处理。这通常比硬判决的Viterbi路径效果更好。 % 软收缩系数 * P(状态大系数) weights gamma(:, 2); % 属于状态2大系数的概率 processed_coeffs_flat_soft obs_sequence .* weights; % 同样需要将processed_coeffs_flat_soft映射回C_recon。这里我们假设能完美映射。 % 重构信号 % 假设我们有一个函数能将处理后的系数向量还原成小波分解结构C_processed % C_processed reverse_mapping(processed_coeffs_flat_soft, C, L, J, wavelet_name); % x_denoised waverec(C_processed, L, wavelet_name); % 由于逆向映射较复杂为了完成仿真闭环我们做一个概念性演示 % 我们直接对原始噪声信号用小波硬阈值去噪并与HMM软收缩的思想对比但非同一基准。 % 使用wden函数进行默认阈值去噪 [x_wden, ~, ~] wden(x_noisy, rigrsure, s, mln, J, wavelet_name); figure; subplot(2,2,1); plot(t, x_clean); title(原始干净信号); grid on; ylim([-1 4]); subplot(2,2,2); plot(t, x_noisy); title(加噪信号); grid on; ylim([-1 4]); subplot(2,2,3); plot(t, x_wden); title([小波默认阈值去噪 (, wavelet_name, )]); grid on; ylim([-1 4]); % 计算并显示一些指标 mse_wden mean((x_wden - x_clean).^2); psnr_wden 10*log10(max(x_clean.^2)/mse_wden); text(0.05, 0.9, sprintf(MSE: %.4f\nPSNR: %.2f dB, mse_wden, psnr_wden), Units, normalized, FontSize, 9); subplot(2,2,4); % 这里本应绘制HMM去噪结果。我们用文本说明替代。 plot(t, x_noisy, Color, [0.8 0.8 0.8]); hold on; title(HMM状态概率去噪 (概念示意)); xlabel(时间); ylabel(幅度); grid on; ylim([-1 4]); % 绘制状态概率作为背景 imagesc([t(1), t(end)], [-1, 4], repmat(weights, 2, 1)); % 需要调整weights的维度以匹配时间 colormap(flipud(gray)); colorbar; clim([0 1]); alpha(0.3); plot(t, x_clean, r-, LineWidth, 1.5, DisplayName, 干净信号); legend(Location, best); text(0.05, 0.9, {灰色背景系数属于大信号状态的概率,红色曲线期望的去噪结果干净信号}, Units, normalized, FontSize, 9);4. 仿真结果分析、常见陷阱与进阶思考运行完上述代码你会得到几个关键的输出EM算法的收敛曲线、估计出的HMM参数A, μ, σ²、每个系数的状态概率γ、以及去噪效果的对比图。4.1 结果解读与模型评估收敛曲线正常情况下对数似然值应该随着迭代单调递增EM算法的性质并逐渐趋于平缓。如果曲线震荡或下降说明代码有bug通常是缩放因子处理不当或概率计算下溢。估计参数转移矩阵A你应该会看到A(1,1)小-小和A(2,2)大-大的概率相对较高这符合信号中“平滑区域”和“边缘区域”各自具有连续性的先验知识。A(1,2)和A(2,1)的概率较低表示状态切换不频繁。高斯参数状态1小系数的估计方差σ²₁应该接近添加的噪声方差0.25均值μ₁接近0。状态2大系数的方差σ²₂会大得多均值μ₂可能偏离0反映了信号突变的幅度。状态概率图在最后的概念图中灰色深浅代表了每个小波系数属于“大信号”状态的概率。你应该能看到在原始干净信号的阶跃和脉冲位置对应的时间区域概率值较高颜色较深而在平坦的噪声区域概率值较低颜色较浅。这直观地展示了WD-HMM如何利用上下文信息跨尺度依赖性来区分信号和噪声比单纯依靠系数幅值的小波阈值法更精准。4.2 仿真中踩过的坑与核心注意事项数值下溢是头号敌人前向概率α_t(i)是许多小于1的数的连乘即使对于中等长度的序列也会迅速变成0。必须实现缩放Scaling。我给出的代码中scale_factor和对应的α、β处理是正确实现的关键。忘记缩放你的EM算法会在前几次迭代后就因为概率全为0而得到NaN参数。观测概率密度B的计算计算高斯概率密度时指数部分(o-μ)²/σ²可能很大导致exp()结果下溢为0。代码中B(B eps) eps;这行是一种简单的防护。更稳健的做法是计算对数概率log B并在前向-后向算法中全程使用对数域计算Log-Sum-Exp技巧这能提供更好的数值稳定性但实现更复杂。对于教学仿真加一个eps下限通常够用。参数初始化至关重要EM算法只能找到局部最优解。糟糕的初始化可能导致收敛到无意义的解例如两个状态的高斯分布参数几乎相同。我的初始化策略基于系数统计粗略划分在简单案例中有效。对于更复杂的信号可以考虑使用K-Means对系数进行聚类用聚类结果初始化μ和σ²用聚类标签的转移频率初始化A。树结构索引映射本次仿真为了突出EM算法核心简化了系数组织使用了从粗到细的简单拼接。这严重削弱了WD-HMM的核心优势——跨尺度依赖性。在真正的WD-HMM中转移概率A是定义在父子节点之间的。你需要构建一个数据结构明确记录每个系数节点的父节点、子节点索引。E步和M步中的求和对t求和需要按照树的拓扑顺序进行。这是将仿真从“玩具”升级到“实用”最关键、也是最繁琐的一步。方差估计的稳定性在M步更新方差时分母是Σ γ_t(k)这是一个介于0和T之间的数。如果某个状态k的后验概率γ_t(k)在所有时刻都很小那么Σ γ_t(k)可能接近0导致方差估计σ_k²异常大。代码中加入max(..., eps)是一种保护。更好的做法是设置一个方差下限比如max(sigma2_new(k), (0.01*coeff_std)^2)。4.3 从仿真到应用模型扩展与优化方向这个基础的WD-HMM仿真框架可以沿多个方向扩展连续密度HMM目前每个状态对应一个单高斯。可以扩展为高斯混合模型GMM即b_k(o) Σ_m c_{km} N(o; μ_{km}, σ_{km}^2)能更精细地刻画一个状态内系数的复杂分布。EM算法需要引入另一个隐变量混合分量指示器更新公式会稍复杂但原理相通。上下文相关建模除了跨尺度依赖还可以建模同尺度内相邻系数兄弟节点的依赖关系这需要更复杂的图模型如马尔科夫随机场估计方法也会变为近似推理如吉布斯采样、置信传播。应用于图像处理将一维信号扩展到二维图像。小波变换变为二维小波树变为四叉树。HMM的状态可以设计为更多种类如水平边缘、垂直边缘、平滑区域等。计算量会显著增加但去噪、分割效果通常比一维更惊艳。与深度学习结合HMM的参数特别是转移矩阵A可以设计成由神经网络预测形成可训练的判别式模型。或者用神经网络来学习观测概率B替代手工设计的GMM。通过这个从理论推导到MATLAB逐行实现的仿真过程我希望你不仅掌握了WD-HMM和EM算法的代码实现更能理解其背后的统计建模思想利用数据的结构先验跨尺度马尔科夫性和概率框架隐变量模型通过迭代优化EM来“学习”数据的生成规律最终实现更智能的信号分析与处理。这比任何黑箱调用现成工具箱的收获都要大得多。下次当你面对一段看似杂乱无章的信号时或许可以想想是不是有一串隐藏的状态链在背后操纵着一切而EM算法正是照亮这条隐藏链条的那束光。