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

news2024/11/25 6:48:01

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

上一篇文章我们讲到了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绘图程序)

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

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

【大厂面试演练】知道ZooKeeper有什么应用场景吗

面试官:咳咳咳,看你简历写了精通ZooKeeper,那我就随便考考你吧 面试官:不用慌尽管说,错了也没关系😊。。。 每日分享大厂面试演练,感兴趣就关注我吧❤️ 面试官:知道ZooKeeper有什么…

Docker安装步骤笔记

一、环境准备 VM网络配置 打开VMware软件 --编辑 --虚拟网络编辑器 二、VM创建虚拟机 三、安装rhel8.9操作系统 1、rhel8.9 镜像下载 第一步:进入redhat官网进行注册第二步:下载rhel8.9镜像文件 https://access.redhat.com/downloads/content/rhel …

Pytorch搭建AlexNet 预测实现

1.导包 import torch import matplotlib.pyplot as plt import json from model import AlexNet from PIL import Image from torchvision import transforms 2.数据预处理 data_transform transforms.Compose([transforms.Resize((224, 224)), # 将图片重新裁剪transform…

JDBC连接Mysql(executeQuely)3/13

resultset-->executeQuery import java.sql.Connection; import java.sql.DriverManager; import java.sql.ResultSet; import java.sql.Statement;public class Demo3 {public static void main(String[] args) throws Exception {//1.注册驱动Class.forName("com.mys…

【Java,Redis】Redis 数据库存取字符串数据以及类数据

1、 字符串存取数据 Resource private StringRedisTemplate stringRedisTemplate;//从Redis中获取string字符串 stringRedisTemplate.opsForValue().get("cache:shop:"id); //Json -> class Shop shop JSONUtil.toBean(ShopJson,Shop.class); //字符串写入redis…

C# Stopwatch计算代码运行时间

文章目录 前言一、计算范围时间1、起始位置2、结束位置3、获取时间封装成对象(1)、完整代码(2)、使用示例 二、计算检查点时间1、初始化2、检查点封装成对象(1)、完整代码(2)、使用示…

SQL Server错误:15404

执行维护计划失败,提示SQL Server Error 15404 无法获取有关... 异常如下图: 原因:数据库用户名与计算机名称不一致 解决办法:1.重名称数据库用户名 将前缀改成计算机名 2.重启SQL Server代理

C++Qt学习——不用UI文件编程

在创建文件的时候不要选中Generate form这块 创建的文件如下图所示,比起之前的没有了form这一快 1、在mainwindow.h里面声明按钮对象 2、在mainwindow.cpp里实例化按钮 2.1、方法一 pushButton new QPushButton();pushButton->show(); 但是发现显示是分离的 2…

保研复习数据结构记(8)--排序

一.内部排序 1.概念 什么是排序?是将一个任意排列的记录或者是数据元素,排列成按关键字有序的序列什么是排序方法是稳定的?按照关键字排序的kikj,在排序之后,两个关键字相等的记录的顺序与排序之前相同,若…

ubuntu2004桌面系统英伟达显卡驱动安装方法

#如何查看显卡型号 lspci | grep -i vga#----output------ 01:00.0 VGA compatible controller: NVIDIA Corporation Device 1f06 (rev a1)根据 Device 后的 值 进入网站查询 pci-ids.ucw.cz/mods/PC/10de?actionhelp?helppci #根据显卡型号,下载对应系统的驱动…

Linux 硬件时间(RTC time)、系统时间(UTC时间、Universal time)、本地时间(Local time)、时区(Time zone)与夏令时(DST)解析

文章目录 理解时间:硬件时间、系统时间(UTC时间)、本地时间、时区与夏令时1. 硬件时间(RTC time)1.1 硬件时间简介1.2 如何使用硬件时间 2. 系统时间(UTC时间)(Universal time&#…

TSINGSEE青犀煤矿矿井视频监控与汇聚融合管理视频监管平台建设方案

一、背景需求 随着我国经济的飞速发展,煤炭作为我国的主要能源之一,其开采和利用的重要性不言而喻。然而,煤矿事故频发,不仅造成了巨大的人员伤亡和财产损失,也对社会产生了深远的负面影响。视频监控系统作为实现煤矿智…

普发Pfeiffer OmniStar/ThermoStar GSD300/GSS300内部电路图装配安装3D图原理图电路板电路图详情内容看图片目录

普发Pfeiffer OmniStar/ThermoStar GSD300/GSS300内部电路图装配安装3D图原理图电路板电路图详情内容看图片目录

【考研】高等数学总结

文章目录 第一章 极限 函数 连续1.1 极限存在准则及两个重要极限1.1.1 夹逼定理1.1.1.1 数列夹逼定理1.1.1.2函数夹逼定理 1.1.2 两个重要极限1.1.2.1 极限公式11.1.2.1.1 证明1.1.2.1.2 数列的单调有界收敛准则1.1.2.1.2.1 二项式定理1.1.2.1.2.2 证明 1.1.2.2 极限公式21.1.2…

在Linux/Ubuntu/Debian中设置字体

下载字体。 下载你喜欢的字体,双击并安装。 之后更新字体缓存: fc-cache -f -v安装 GNOME 调整。 GNOME Tweaks 是一个工具,允许你自定义 GNOME 桌面环境的各个方面,包括字体。 如果你还没有安装 GNOME Tweaks: …

pytorch之诗词生成--2

先上代码: # -*- coding: utf-8 -*- # File : dataset.py # Author : AaronJny # Time : 2019/12/30 # Desc : 构建数据集 from collections import Counter import math import numpy as np import tensorflow as tf import settingsclass Tokenizer:""&…

MIT 6.S081---Lab: locks

Memory allocator (moderate) 修改kernel/kalloc.c,修改kmem声明并定义结构体数组: 修改kernel/kalloc.c中的kinit函数,对kmemList进行初始化: 修改kernel/kalloc.c中的kfree函数,获取当前的cpuid并将释放的内存添加到…

互联网架构与通信机制:从边缘到核心的深度解析

✨✨ 欢迎大家来访Srlua的博文(づ ̄3 ̄)づ╭❤~✨✨ 🌟🌟 欢迎各位亲爱的读者,感谢你们抽出宝贵的时间来阅读我的文章。 我是Srlua小谢,在这里我会分享我的知识和经验。&am…

vscode使用npm命令无反应,而终端可以的解决办法

如若你遇到这种情况 使用命令 get-command npm 去下面这个路径把它删掉就可以了

MyBatis拦截器四种类型和自定义拦截器的使用流程

文章目录 MyBatis拦截器四种类型和自定义拦截器的使用流程一、MyBatis拦截器四种类型的详细解释:1. **ParameterHandler 拦截器**:2. **ResultSetHandler 拦截器**:3. **StatementHandler 拦截器**:4. **Interceptor Chain 拦截器…