【滤波专题-第8篇】ICA降噪方法——类EMD联合ICA降噪及MATLAB代码实现（以VMD-ICA为例）

今天来介绍一种效果颇为不错的降噪方法。（针对高频白噪声）

上一篇文章我们讲到了FastICA方法。在现实世界的许多情况下，噪声往往接近高斯分布，而有用的信号（如语音、图像特征等）往往表现出非高斯的特性。FastICA通过最大化输出信号的非高斯性来恢复这些有用的信号，从而有效地从噪声中分离出信号。

将类经验模态分解（EMD）方法与FastICA结合使用，可以创建一种巧妙的信号处理策略，这种结合利用了两种方法的互补优势：通过EMD分解得到的IMFs可以简化信号的结构，但单独使用EMD可能不足以有效分离信号中的噪声和有用信息，自适应得到的imf数量也可能存在冗余；FastICA则能够从imf信号中进一步提取独立成分，更有效地识别噪声成分。

下面将详细解释这种结合的算法流程、优势以及MATLAB代码实现。

一、算法流程

我们直接开门见山，通过案例来讲解算法流程。

首先生成一段待降噪信号。

rng(123456)
t = (0:0.001:(1-0.001))';
x1 = 0.6*sin(15*pi*t+pi/5);
x2 = cos(60*pi*t+sin(10*pi*t));
x3 = (1+0.3*cos(10*pi*t)).*sin(200*pi*t);
x4 = wgn(1000,1,-10);

fs = 1000; %采样频率
ps = x1+x2+x3;  %未加噪声的纯净信号
s = ps+x4;
figure('color','w')
subplot(2,1,1);plot(t,ps,'k');title('无噪声的仿真信号')
subplot(2,1,2);plot(t,s,'k');title('加入噪声的仿真信号')

这段信号加入噪声前后的图像如下，我们的任务就是将下边的含噪信号尽可能回复成上图中的无噪声信号。

1.类EMD分解

关注我的专栏的老读者们都知道，我将EMD、EEMD、CEEMD、CEEMDAN、ICEEMDAN、VMD等一系列模态方法统称为“类EMD”方法，在本篇文章中，该步骤也可以是上述分解方法中的任意一种，毕竟这一步的目的是将复杂信号分解为一定数量的模态分量。

此案例中我们使用VMD分解算法。

alpha=1000;  % alpha   - 惩罚因子
tol=1e-6;    % tol     - 收敛容差，是优化的停止准则之一，可以取 1e-6~5e-6
K=5;         % K       - 指定分解模态数
imf = pVMD(s,fs, alpha, K, tol);

上边代码中使用了pVMD函数，这是笔者封装的一个易用的VMD画图函数，其介绍可以看之前的文章。分解结果如下：

也可以在调用pVMDandFFT函数绘制imf分量与频谱对比图。

imf = pVMDandFFT(s,fs, alpha, K, tol); %函数见此处：https://zhuanlan.zhihu.com/p/396775790

当然使用MATLAB内置的vmd函数也是可以的，但是要注意输入FastICA的imf分量的行列方向，需要保证imf是每一行为一个分量。

2.FastICA盲源分离

以类EMD分解得到的imf分量作为输入，进行ICA混解。

在这里一定要算出混合矩阵A（忘记概念的同学可以看上一篇）。

%% 3.ICA混解并绘制混解后的图像
numOfIC = 0;  % 需要提取的独立成分数目，如果不指定数目，则输入0
g = 'pow3';   % 使用的非线性函数类型，可选'pow3', 'tanh', 'gauss'
[icasig, A, W] = pFastICA(imf, numOfIC, g);

此时，我们将会得到如下分解结果：

从FastICA盲源分离结果可以推断独立成分5是噪声分量，但是在自动化实现降噪的代码中，我们自然希望对噪声分量的判断是不需要人为介入的。此时我们引入一个功率谱熵阈值判断方法来实现。

3.功率谱熵阈值判断

在之前的文章（Mr.看海：【熵与特征提取】基于“信息熵”的特征指标及其MATLAB代码实现（功率谱熵、奇异谱熵、能量熵））中讲到，信息熵越大，代表不确定性越大，信号中包含的信息量越少，信号越无秩序/越接近于白噪声，则信息熵越大。

由于盲源分离本身的特性，分离出来的独立成分中，最多包含一个高斯分量，也就是噪声分量。

所以我们只需要判断几个独立成分的功率谱熵中，是否有数值相较于其他分量大得多的成分，就是噪声分量。

也就是找到功率谱熵中的离群值。

这里我们阈值的计算标准是：阈值=功率谱熵均值+功率谱熵的标准差

我们将大于阈值的独立成分直接赋值为0。

计算结果如下图，可以看到噪声分量被识别了出来。

4.重构滤波信号

不知道大家是否还记得FastICA方法的重要性质之一：输出信号幅度的不确定性。FastICA分解得到的独立成分是不能像类EMD分解得到的imf那样直接相加进行重构的。

但是FastICA也有自己的重构方法，就是使用上边提到的混合矩阵A。重构公式为：

X=A∗icasig

上式中的icasig就是经过处理过后的独立成分，对于大于阈值的成分，我们已经将其幅值全部置0；X就是重构结果，也就是该方法中的滤波结果。

此时得到的滤波结果为：

可以看出，滤波效果还是颇为不错的，经计算此时的信噪比为15.9dB。

二、算法流程图

将上述算法梳理成流程图如下，当然vmd分解这步是可以替换成其他分解算法的，另外功率谱熵也可以换成其他熵值指标。

三、MATLAB实现

为了方便大家使用，按照惯例笔者对上述流程进行了封装。封装的流程就是上述流程图中红框内的部分，之所以没有将vmd也一并封装进去，是为了大家方便替换其他的模态分解算法。

封装后的函数只需要将分解好的imf变量导入，即可一行代码实现滤波：

reSig = filEMDsICA(imf);  % 调用filEMDsICA函数实现滤波

经实测，这个方法的滤波效果还是很不错的。

我为大家提供了完整的演示案例，可以一键实现下述图像和滤波结果的输出：

要你做的，基本就只需替换数据啦！

需要上边这个函数文件以及测试代码的同学，可以在公众号 khscience（看海的城堡）中回复“ICA滤波”获取。

扩展阅读

Mr.看海：【滤波专题-第1篇】数字滤波器15分钟入门！——这可能是最简单的FIR有限冲激响应滤波讲解

Mr.看海：【滤波专题-第2篇】数字滤波器15分钟入门！——这可能是最简单的IIR无限冲激响应滤波讲解

Mr.看海：【滤波专题-第3篇】IIR无限冲激响应和FIR有限冲激响应数字滤波器有什么区别？

Mr.看海：【滤波专题-第4篇】滤波器滤波效果的评价指标（信噪比SNR、均方误差MSE、波形相似参数NCC）

Mr.看海：【滤波专题-第5篇】FIR、IIR滤波器设计及MATLAB实现

Mr.看海：【滤波专题-第6篇】小波阈值去噪方法看这一篇就明白了~（附MATLAB实现）

Mr.看海：【滤波专题-第7篇】“类EMD”算法分解后要怎样使用（3）——EMD降噪方法及MATLAB代码实现

Mr.看海：【盲源分离】快速理解FastICA算法（附MATLAB绘图程序）