Matlab小波去噪——基于wden函数的去噪分析

news2024/12/28 20:06:20

文章目录

  • 一、问题描述
  • 二、代码
    • 问题1:原始信号加6分贝高斯白噪声
    • 问题2:确定合适的小波基函数
    • 问题3:确定最合适的阈值计算估计方法
    • 问题4:确定合适的分解层数
    • 问题5:实际信号去噪
    • 问题6:对比
  • 三、演示视频
  • 最后


一、问题描述

1.利用MATLAB绘制原始信号,对其加6分贝高斯白噪声;
2.以Minimaxi阈值法,软阈值函数,3层分解层数,分别用dbN和symN小波对加噪信号去噪,获得分解图和去噪后的图,并用信噪比和均方根误差作为评判标准,确定合适的小波基函数;
3.用第2步确定的小波基函数,软阈值函数,分解层数为3层,对无偏估计阈值(RigrSure)、固定式阈值(Sqtwolog)、启发式阈值(HeurSure)和极大极小阈值(Minimaxi)四种分别去噪,获得去噪后的图,并用信噪比和均方根误差作为评判标准,确定最合适的阈值计算估计方法;
4.用第2步确定的小波基函数,第3步确定的阈值计算估计准则,分别用分解层数为1,2,3,4,5,6对加噪信号进行去噪,获得去噪后得到图,并用信噪比和均方根误差作为评判标准;
5.用实际的信号加6分贝噪声对前面确定的小波基函数,阈值计算方法以及分解层数用小波阈值进行去噪,并求信噪比和均方根误差。
6、确定好小波基函数、阈值函数和分解层数后,分别模拟加入不同量的噪声与4阶巴特沃斯低通滤波器滤波对比

二、代码

问题1:原始信号加6分贝高斯白噪声

代码如下(示例):

clear
clc
close all
%% MATLAB绘制原始信号
load('data.mat'); %私聊发数据
data=data;
%%6分贝高斯白噪声
SNR=6; %6dB
noise=0.2*randn(size(data))*std(data)/db2mag(SNR);
s=data+noise;
figure;
subplot(211)
plot(data);ylabel('P/MPa');title('原始信号')
subplot(212)
plot(s);ylabel('P/MPa');title('加6dB高斯白噪声')

问题2:确定合适的小波基函数

代码如下(示例):

clear
clc
close all
%% MATLAB绘制原始信号
load('data.mat');
data=data;
%%6分贝高斯白噪声
SNR=6; %6dB
noise=0.2*randn(size(data))*std(data)/db2mag(SNR);
s=data+noise;
%% Minimaxi阈值法,软阈值函数,3层分解层数,db5去噪
wname=strvcat('sym4','sym5','db4','db5');
for i=1:4
    [C,L] = wavedec(s,3,wname(i,:));  %进行3层小波包分解
    s1=wden(s,'minimaxi','s','mln',3,wname(i,:)); %Minimaxi、软阈值,3层,db5
    figure;
    subplot(311)
    plot(data);xlabel('t/ms');ylabel('P/MPa');title('原始信号')
    subplot(312)
    plot(s);xlabel('t/ms');ylabel('P/MPa');title('加6dB高斯白噪声')
    subplot(313)
    plot(s1);xlabel('t/ms');ylabel('P/MPa');title(['Minimaxi-软阈值-3层-',wname(i,:)])
    figure
    subplot(511)
    plot(data,'r');ylabel('s');title([wname(i,:),'小波分解图'])
    set(gca,'ytick',[]) 
    set(gca,'xtick',[]) 
    subplot(512)
    plot(C(1:L(2)),'b');ylabel('a3')
    set(gca,'ytick',[]) 
    set(gca,'xtick',[]) 
    subplot(513)
    plot(C(L(2):L(3)));ylabel('d3')
    set(gca,'ytick',[]) 
    set(gca,'xtick',[]) 
    subplot(514)
    plot(C(L(3):L(4)));ylabel('d2')
    set(gca,'ytick',[]) 
    set(gca,'xtick',[]) 
    subplot(515)
    plot(C(L(4):L(5)));ylabel('d1')
    SNR_s1(i)=snr(data,s1-data);RMSE_s1(i)=sqrt(mse(data-s1));   
    SNR_s11(i)=snr(s,s1-s);
    disp(['Minimaxi-软阈值-3层-',wname(i,:),':信噪比=',num2str(SNR_s1(i)),'dB,均方根误差=',num2str(RMSE_s1(i))])
    disp(['加噪后信噪比=',num2str(SNR_s11(i)),'dB'])
    disp('-----------------------------------------------------------')
end
%% 根据SNR选取较好的小波基函数
[m,index]=max(SNR_s1);
disp(['最合适的阈值计算估计方法为:',wname(index,:)])
disp('-----------------------------------------------------------')

问题3:确定最合适的阈值计算估计方法

代码如下(示例):

clear
clc
close all
%% MATLAB绘制原始信号
load('data.mat');
data=data;
%%6分贝高斯白噪声
SNR=6; %6dB
noise=0.2*randn(size(data))*std(data)/db2mag(SNR);
s=data+noise;
%% main2已经确定最合适的小波基函数
wname='sym5';
%% 无偏估计阈值(RigrSure)、固定式阈值(Sqtwolog)、启发式阈值(HeurSure)和极大极小阈值(Minimaxi)
TPTR=['rigrsure';'sqtwolog';'heursure';'minimaxi'];
for i=1:4
    s3=wden(s,TPTR(i,:),'s','mln',3,wname); %依次进行滤波
    figure
    subplot(311)
    plot(data);xlabel('t/ms');ylabel('P/MPa');title('原始信号')
    subplot(312)
    plot(s);xlabel('t/ms');ylabel('P/MPa');title('加6dB高斯白噪声')
    subplot(313)
    plot(s3);xlabel('t/ms');ylabel('P/MPa');title(['采用',TPTR(i,:),'进行滤波'])
    snr_s3(i)=snr(data,s3-data);RMSE_s3(i)=sqrt(mse(data-s3));
    snr_s33(i)=snr(s,s3-s);
    disp([TPTR(i,:),'-软阈值-3层-',wname,':信噪比=',num2str(snr_s3(i)),'dB,均方根误差=',num2str(RMSE_s3(i))])
    disp(['加噪后信噪比=',num2str(snr_s33(i)),'dB'])
    disp('-----------------------------------------------------------')
end
%% 根据SNR选取较好的阈值计算估计方法
[m,index]=max(snr_s3);
disp(['最合适的阈值计算估计方法为:',TPTR(index,:)])
disp('-----------------------------------------------------------')

问题4:确定合适的分解层数

代码如下(示例):

clear
clc
close all
%% MATLAB绘制原始信号
load('data.mat');
data=data;
%%6分贝高斯白噪声
SNR=6; %6dB
noise=0.2*randn(size(data))*std(data)/db2mag(SNR);
s=data+noise;
%% main2和main3确定的小波基函数和阈值计算估计方法
wname='sym5';
TPTR='sqtwolog';
%% 分解层数为123456
for i=1:6
    s4=wden(s,TPTR,'s','mln',i,wname); %依次进行滤波
    figure
    subplot(311)
    plot(data);xlabel('t/ms');ylabel('P/MPa');title('原始信号')
    subplot(312)
    plot(s);xlabel('t/ms');ylabel('P/MPa');title('加6dB高斯白噪声')
    subplot(313)
    plot(s4);xlabel('t/ms');ylabel('P/MPa');title(['分解层数=',num2str(i)])
    snr_s4(i)=snr(data,s4-data);RMSE_s4(i)=sqrt(mse(data-s4));
    snr_s44(i)=snr(s,s4-s);
    disp([TPTR,'-软阈值-',num2str(i),'层-',wname,':信噪比=',num2str(snr_s4(i)),'dB,均方根误差=',num2str(RMSE_s4(i))])
    disp(['加噪后信噪比=',num2str(snr_s44(i)),'dB'])
    disp('-----------------------------------------------------------')
end
%% 根据SNR选取较好的分解层数
[m,index]=max(snr_s4);
disp(['最合适的分解层数为:',num2str(index)])
disp('-----------------------------------------------------------')

问题5:实际信号去噪

代码如下(示例):

clear
clc
close all
%% 读取实际的信号
data=xlsread('14#c1.csv');
data=data(:,2);
%%6分贝高斯白噪声
SNR=6; %6dB
noise=0.2*randn(size(data))*std(data)/db2mag(SNR);
s=data+noise;
%% 根据(2)(3)(4)确定参数
wname='sym5';
TPTR='sqtwolog';
lev=6;
%% 进行滤波
s5=wden(s,TPTR,'s','mln',lev,wname); %进行滤波
%% 绘制
figure;
subplot(311)
plot(data);xlabel('t/ms');ylabel('P/MPa');title('实际信号')
subplot(312)
plot(s);xlabel('t/ms');ylabel('P/MPa');title('加6dB高斯白噪声')
subplot(313)
plot(s5);xlabel('t/ms');ylabel('P/MPa');title('信号去噪')
snr_s55=snr(s,s5-s);
snr_s5=snr(data,s5-data);RMSE_s5=sqrt(mse(data-s5));
disp([TPTR,'-软阈值-',num2str(lev),'层-',wname,':信噪比=',num2str(snr_s5),'dB,均方根误差=',num2str(RMSE_s5)])
disp(['加噪后信噪比=',num2str(snr_s55),'dB'])

问题6:对比

代码如下(示例):

clear
clc
close all
%% 读取实际的信号
data=xlsread('14#c1.csv');
data=data(:,2);
fs=125000; 
%%
wname='sym5';
TPTR='sqtwolog';
lev=6;
%% 设计4阶巴特沃斯低通滤波器
fc=10000;
n=4;  %阶数
[b,a]=butter(n,fc/(fs/2), 'low');
%%1-16分贝高斯白噪声
for SNR=1:16 %6dB
    noise=0.2*randn(size(data))*std(data)/db2mag(SNR);
    s=data+noise;
    s1=filter(b,a,s);  %filter既能进行IIR滤波又能进行FIR滤波
    s2=wden(s,TPTR,'s','mln',lev,wname); %进行滤波
    snr_s1(SNR)=snr(data,s1-data);RMSE_s1(SNR)=sqrt(mse(data-s1));
    snr_s2(SNR)=snr(data,s2-data);RMSE_s2(SNR)=sqrt(mse(data-s2));
end
figure;
plot(1:16,snr_s1,'o-r');
hold on
plot(1:16,snr_s2,'*-b');
xlabel('高斯白噪声dB');ylabel('SNR')
legend('FIR滤波','小波滤波')
title('信噪比曲线')
%% 
for SNR=2:2:10 %6dB
    noise=0.2*randn(size(data))*std(data)/db2mag(SNR);
    s=data+noise;
    s1=filter(b,a,s);  %filter既能进行IIR滤波又能进行FIR滤波
    s2=wden(s,TPTR,'s','mln',lev,wname); %进行滤波
    figure;
    subplot(411)
    plot(data);xlabel('t/ms');ylabel('P/MPa');title('实际信号')
    subplot(412)
    plot(s);xlabel('t/ms');ylabel('P/MPa');title('加16dB高斯白噪声')
    subplot(413)
    plot(s1);xlabel('t/ms');ylabel('P/MPa');title('FIR信号去噪')
    subplot(414)
    plot(s1);xlabel('t/ms');ylabel('P/MPa');title('小波信号去噪')
    suptitle(['噪声大小=',num2str(SNR),'dB'])
end


三、演示视频

基于wden函数的去噪演示

最后

不定期发布一些matlab设计内容,敬请期待。包括但不限于如下内容:信号处理、通信仿真、gui设计、matlab appdesigner,simulink仿真。有任何有关MATLAB的问题可🐧咨询
在这里插入图片描述

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

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

相关文章

团队死气沉沉?10种玩法激活你的项目团队拥有超强凝聚力

作为项目经理和PMO,以及管理者最头疼的是团队的氛围和凝聚力,经常会发现团队死气沉沉,默不作声,你想尽办法也不能激活团队,也很难凝聚团队。这样的项目团队你很难带领大家去打胜仗,攻克堡垒。但是如何才能避…

Python|贪心|数组|二分查找|贪心|数学|树|二叉搜索树|在排序数组中查找元素的第一个和最后一个位置|计数质数 |将有序数组转换为二叉搜索树

1、在排序数组中查找元素的第一个和最后一个位置(数组,二分查找) 给定一个按照升序排列的整数数组 nums,和一个目标值 target。找出给定目标值在数组中的开始位置和结束位置。 如果数组中不存在目标值 target,返回 […

第十四届蓝桥杯三月真题刷题训练——第 2 天

目录 题目1:奇数倍数 代码: 题目2:求值 代码: 题目3:求和 代码: 题目4:数位排序 代码: 题目1:奇数倍数 题目描述 本题为填空题,只需要算出结果后,在代码中使用输出语句将所填结果输出即…

收银系统的设计与实现

技术:Java、JSP等摘要:随着销售行业竞争的日益激烈,收银系统的引入显得极其重要。收银系统不但可以提高商品存储管理的工作效率,而且可以有效减少盲目采购、降低采购成本、合理控制库存、减少资金占用并提高市场灵敏度&#xff0c…

Java虚拟机的运行时数据区-go语言实现

Java虚拟机的运行时数据区 Java虚拟机把存放各式各样数据的内存区域叫作运行时数据区。运行时数据区分成两类: 一类时多线程共享的,一类时线程私有的。多线程共享的数据在Java虚拟机启动时创建好,在Java虚拟机退出时销毁。线程私有的运行时…

序列号和反序列化--java--Serializable接口--json序列化普通使用

序列化和反序列化序列化和反序列化作用为什么需要用途Serializable使用serialVersionUID不设置的后果什么时候修改Externalizable序列化的顺序json序列化序列化和反序列化 序列化:把对象转换为字节序列的过程称为对象的序列化。 反序列化:把字节序列恢复为对象的过…

【Go语言学习】安装与配置

文章目录前言一、Go语言学习站二、安装与配置1.安装2.环境变量配置3.Gland编辑器安装与配置Hello, World!总结前言 Go语言特性 Go,又称为 Golang,是一门开源的编程语言,由 Google 开发。Go 语言的设计目标是提供一种简单、快速、高效、安全…

在MySQL中使用不等于符号还能走索引吗?

一般情况下,我们会在一个索引上较多的使用等值查询或者范围查询,此时索引大多可以帮助我们极快的查询出我们需要的数据。 那当我们在where条件中对索引列使用!查询,索引还能发挥他的作用吗? 以此SQL为例: select * …

农产品销售系统的设计与实现

技术:Java、JSP等摘要:这篇文章主要描述的是农产品蔬菜在线销售系统的设计与实现。主要应用关于JSP网站开发技术,并联系到网站所处理的数据的结构特点和所学到的知识,应用的主要是Mysql数据库系统。系统实现了网站的基本功能&…

计算机组成原理|第一章(笔记)

目录第一章 计算机系统概论1.1 计算机系统简介1.1.1 计算机的软硬件概念1.1.2 计算机系统的层次结构1.1.3 计算机组成和计算机体系结构1.2 计算机的基本组成1.2.1 冯 诺伊曼计算机的特点1.2.2 计算机的硬件框图1.2.3 计算机的工作过程1.3 计算机硬件的主要技术指标1.3.1 机器字…

kaggle数据集下载当中所遇到的问题

kaggle数据集下载当中所遇到的问题报错分析pip install kagglethe SSL module is not available解决方法pip的版本升级解决办法下载kaggle包kaggle数据集下载问题解决参考内容报错分析 今天在尝试使用pip install kaggle的方法去下载我需要的数据集的时候遇到了一些报错的问题…

二分查找与判定树

二分查找的算法思想二分查找也称“折半查找”,要求查找表为采用顺序存储结构的有序表。本例一律采用升序排列。二分查找每一次都会比较给定值与序列[low,high]的中间元素,该元素的下标为mid (lowhigh)/2,若两者相等,则返回元素的下标为mid;如…

Django的DRF从入门到精通

第一讲:建立纯净版Django项目 ① 创建Django项目 ② 创建app一个 python manage.py startapp APP名字 ③ 在settings里配置rest_framework,把不需要的全部注释掉 INSTALLED_APPS = [# django.contrib.admin,# django.contrib.auth,# django.contrib.contenttypes,# djang

centos7 安装 MySQL5.7

1.下载MySQL官方的 Yum Repository wget -i -c http://dev.mysql.com/get/mysql57-community-release-el7-10.noarch.rpm2.安装 Yum Repository yum -y install mysql57-community-release-el7-10.noarch.rpm3 使用 yum 安装 MySQL yum -y install mysql-community-server若…

推荐系统1--Deepfm学习笔记

目录 1 keras实现Deepfm demo 2 deepctr模版 3 其他实现方式 ctr_Kera 模型 数据集 预处理 执行步骤 4何为focal loss 参考 1 keras实现Deepfm 假设我们有两种 field 的特征,连续型和离散型,连续型 field 一般不做处理沿用原值,离散型一…

Promise学习基础学习 promise封装fs模块、AJAX请求

Promise 是什么? 抽象表达: 1、Promise 是一门新的技术(ES6规范) 2、Promise 是JS中进行异步编程的新解决方案 备注:旧方案是单纯使用回调函数 具体表达: 1、从语法上来说:Promise 是一个构造…

QML Loader(加载程序)

Loader加载器用于动态加载 QML 组件。加载程序可以加载 QML 文件(使用 source 属性)或组件对象(使用 sourceComponent 属性) 常用属性: active 活动asynchronous异步,默认为falseitem项目progress 进度so…

package.json中 版本号详解

1. 版本号简介 软件版本号有四部分组成: 第一部分:主版本号,当进行不兼容的 API 更改时,则升级主版本;第二部分:次版本号,当以向后兼容的方式添加功能时,则升级次版本;…

FPGA实现SDI视频编解码 SDI接收发送,提供2套工程源码和技术支持

目录1、前言2、设计思路和框架SDI接收SDI缓存写方式处理SDI缓存读方式处理SDI缓存的目的SDI发送3、工程1详解4、工程2详解5、上板调试验证并演示6、福利:工程代码的获取1、前言 FPGA实现SDI视频编解码目前有两种方案: 一是使用专用编解码芯片&#xff0…

【玩转c++】vector讲解和模拟底层实现

本期主题:vector的讲解和模拟实现博客主页:小峰同学分享小编的在Linux中学习到的知识和遇到的问题小编的能力有限,出现错误希望大家不吝赐vector的介绍及使用1.1vector的介绍vector其实就是一个数组的模板 ,存放的数据可以改变而已…