【滤波专题-第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绘图程序)

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/461382.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

unity学习(60)——选择角色界面--MapHandler2-MapHandler.cs

1.新建一个脚本&#xff0c;里面有static变量loadingPlayerList using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Threading.Tasks;namespace Assets.Scripts.Model {internal class LoadData{public static List<Pl…

3D地图在BI大屏中的应用实践

前言 随着商业智能的不断发展&#xff0c;数据可视化已成为一项重要工具&#xff0c;有助于用户更好地理解数据和分析结果。其中&#xff0c;3D地图作为一种可视化工具&#xff0c;已经在BI大屏中得到了广泛地应用。 3D地图通过将地理信息与数据相结合&#xff0c;以更加直观…

【AI】用iOS的ML(机器学习)创建自己的AI App

用iOS的ML(机器学习)创建自己的AI App 目录 用iOS的ML(机器学习)创建自己的AI App机器学习如同迭代过程CoreML 的使用方法?软件要求硬件开始吧!!构建管道:设计和训练网络Keras 转 CoreML将模型集成到 Xcode 中结论推荐超级课程: Docker快速入门到精通Kubernetes入门到…

计算机网络——物理层(数据通信基础知识)

计算机网络——物理层&#xff08;1&#xff09; 物理层的基本概念数据通信的基本知识一些专业术语消息和数据信号码元 传输速率的两种表示方法带宽串行传输和并行传输同步传输和异步传输 信道基带信号调制常用编码方式 我们今天进入物理层的学习&#xff0c;如果还没有了解OSI…

Transformer代码从零解读【Pytorch官方版本】

文章目录 1、Transformer大致有3大应用2、Transformer的整体结构图3、如何处理batch-size句子长度不一致问题4、MultiHeadAttention&#xff08;多头注意力机制&#xff09;5、前馈神经网络6、Encoder中的输入masked7、完整代码补充知识&#xff1a; 1、Transformer大致有3大应…

C++ 入门篇

目录 1、了解C 2、C关键字 2、命名空间 2.1 命名空间的定义 2.2 命名空间的使用 3. C输入与输出 4.缺省参数 4.1 缺省参数的概念 4.2 缺省参数的分类 5. 函数重载 5.1 函数重载的概念 5.2 C中支持函数重载的原理--名字修饰 6. 引用 6.1 引用概念 6.2 引用…

Docker 系列2【docker安装mysql】【开启远程连接】

文章目录 前言开始步骤1.增加mysql挂载目录2.下载镜像2.启动容器具体步骤4.无法连接5.测试连接 总结 前言 本文开始&#xff0c;默认已经安装docker&#xff0c;如果你还没有完成这个步骤&#xff0c;请查看这一篇文章【docker安装与使用】 开始步骤 1.增加mysql挂载目录 m…

网络原理(1)——UDP协议

目录 一、应用层 举个例子&#xff1a;点外卖 约定数据格式简单粗暴的例子 客户端和服务器的交互&#xff1a; 序列化和返序列化 xml、json、protobuffer 1、xml 2、json 3、protobuffer 二、传输层 端口 端口号范围划分 认识知名的端口号 三、UDP协议 端口 U…

宜搭faas服务器报错Network response was not OK

[error] https://api.dingtalk.com/v1.0/yida/forms/instances? fetch error Error: Network response was not OK 不出意外的话肯定是请求代码的某个部分出了问题&#xff1a;其中formInstanceId和updateFormDataJson是业务的内容 我检查过是没问题的。appType和systemToken…

面试经典-MySQL篇

一、MySQL组成 MySQL数据库的连接池&#xff1a;由一个线程来监听一个连接上请求以及读取请求数据&#xff0c;解析出来一条我们发送过去的SQL语句SQL接口&#xff1a;负责处理接收到的SQL语句查询解析器&#xff1a;让MySQL能看懂SQL语句查询优化器&#xff1a;选择最优的查询…

【C#】【SAP2000】读取SAP2000中所有Frame对象在指定工况的温度荷载值到Grasshopper中

if (build true) {// 连接到正在运行的 SAP2000// 使用 COM 接口获取 SAP2000 的 API 对象cOAPI mySapObject (cOAPI)System.Runtime.InteropServices.Marshal.GetActiveObject("CSI.SAP2000.API.SapObject");// 获取 SAP2000 模型对象cSapModel mySapModel mySap…

Vue 项目安装依赖提示core-js版本低的处理办法

core-js2.6.12: core-js<3 is no longer maintained and not recommended for usage due to the number of issues. Please, upgrade your dependencies to the actual version of core-js3. 我是下载一个老的项目&#xff0c;npm install之后提示上面的错误&#xff1b;本…

Linux——ELK日志分析系统

实验环境 虚拟机三台CentOS 7.9&#xff0c; 组件包 elasticsearch-5.5.0.rpm elasticsearch-head.tar.gz node-v8.2.1.tar.gz phantomjs-2.1.1-linux-x86_64.tar.bz2 logstash-5.5.1.rpm kibana-5.5.1-x86_64.rpm 初始…

分享一下自己总结的7万多字java面试笔记和一些面试视频,简历啥的,已大厂上岸

分享一下自己总结的7万多字java面试笔记和一些面试视频&#xff0c;简历啥的&#xff0c;已大厂上岸 自己总结的面试简历资料&#xff1a;https://pan.quark.cn/s/8b602fe53b58 文章目录 SSMspringspring 的优点&#xff1f;IoC和AOP的理解**Bean 的生命周期****列举一些重要…

C++进阶:详解多态(多态、虚函数、抽象类以及虚函数原理详解)

C进阶&#xff1a;详解多态&#xff08;多态、虚函数、抽象类以及虚函数原理详解&#xff09; 结束了继承的介绍&#xff1a;C进阶&#xff1a;详细讲解继承 那紧接着的肯定就是多态啦 文章目录 1.多态的概念2.多态的定义和实现2.1多态的构成条件2.2虚函数2.2.1虚函数的概念2…

算法笔记 连载中。。。

HashMap&#xff08;会根据key值自动排序&#xff09; HashMap<String, Integer> hash new HashMap<>() hash.put(15,18) hash.getOrDefault(ts, -1) //如果ts(key)存在&#xff0c;返回对应的value 否则返回-1 hashMap1.get(words1[i])1会报错&#xff0c;因…

Vue2在一个页面内动态切换菜单显示对应的路由组件

项目的需求是在一个页面内动态获取导航菜单&#xff0c;导航菜单切换的时候显示对应的路由页面&#xff0c;类似于tab切换的形式&#xff0c;切换的导航菜单和页面左侧导航菜单是同一个路由组件&#xff0c;只是放到了一个页面上&#xff0c;显示的个数不同&#xff0c;所有是动…

QT下跨平台库实现及移植经验分享

最近在移植公司一个QT桌面软件到android上&#xff0c;有一些公司自定义的库&#xff0c;用了很多windows的api&#xff0c;移植过程很是曲折&#xff0c;在此有一些感悟分享一下~ 一.自编写跨平台库 1.有时候为了程序给第三方用需要编译一些qt封装库&#xff0c;并可能跨平台…

AI智慧校园电子班牌云平台源码

目录 家长端 学校端 电子围栏 亲情通话 课堂答题 移动化管理模式 统一资源管理平台 模板内容智能更换 家校互联 家长端 多场景通话:上学放学联系、紧急遇险求助联系、日常亲情通话关注孩子人身安全:到校离校情况、进入危险区域预警等。 学校端 课堂秩序管理:提高教…

canvas绘制时,画布上有一个镂空的圆形(即背景可见),然后随着动画的进行,这个圆形的边界逐渐扩大至充满整个屏幕

<canvas id"myCanvas" width"800" height"600"></canvas>在不同宽高比的屏幕上&#xff0c;如果canvas元素没有被强制保持与窗口同样的宽高比&#xff08;例如通过CSS设置其宽度和高度百分比或者响应式布局&#xff09;&#xff0c;…