python小波去噪的方法_小波去噪基本概念
一、前言
在现实生活和工作中,噪声无处不在,在许多领域中,如天文、医学图像和计算机视觉方面收集到的数据常常是含有噪声的。噪声可能来自获取数据的过程,也可能来自环境影响。由于种种原因,总会存在噪声,噪声的存在往往会掩盖信号本身所要表现的信息,所以在实际的信号处理中,常常需要对信号进行预处理,而预处理最主要的一个步骤就是降噪。
小波分析是近年来发展起来的一种新的信号处理工具,这种方法源于傅立叶分析,小波(wavelet),即小区域的波,仅仅在非常有限的一段区间有非零值,而不是像正弦波和余弦波那样无始无终。小波可以沿时间轴前后平移,也可按比例伸展和压缩以获取低频和高频小波,构造好的小波函数可以用于滤波或压缩信号,从而可以提取出已含噪声信号中的有用信号。
二、小波去噪的原理
Donoho提出的小波阀值去噪的基本思想是将信号通过小波变换(采用Mallat算法)后,信号产生的小波系数含有信号的重要信息,将信号经小波分解后小波系数较大,噪声的小波系数较小,并且噪声的小波系数要小于信号的小波系数,通过选取一个合适的阀值,大于阀值的小波系数被认为是有信号产生的,应予以保留,小于阀值的则认为是噪声产生的,置为零从而达到去噪的目的。
从信号学的角度看 ,小波去噪是一个信号滤波的问题。尽管在很大程度上小波去噪可以看成是低通滤波 ,但由于在去噪后 ,还能成功地保留信号特征 ,所以在这一点上又优于传统的低通滤波器。由此可见 ,小波去噪实际上是特征提取和低通滤波的综合 ,其流程图如下所示:
一个含噪的模型可以表示如下:
其中 ,f( k)为有用信号,s(k)为含噪声信号,e(k)为噪声,ε为噪声系数的标准偏差。
假设,e(k)为高斯白噪声,通常情况下有用信号表现为低频部分或是一些比较平稳的信号,而噪声信号则表现为高频的信号,我们对 s(k)信号进行小波分解的时候,则噪声部分通常包含在HL、LH、HH中,如下图所示,只要对HL、LH、HH作相应的小波系数处理,然后对信号进行重构即可以达到消噪的目的。
我们可以看到,小波去噪的原理是比较简单类,类似以往我们常见的低通滤波器的方法,但是由于小波去找保留了特征提取的部分,所以性能上是优于传统的去噪方法的。
三、小波去噪的基本方法
一般来说, 一维信号的降噪过程可以分为 3个步骤
信号的小波分解。选择一个小波并确定一个小波分解的层次N,然后对信号进行N层小波分解计算。
小波分解高频系数的阈值量化。对第1层到第N层的每一层高频系数(三个方向), 选择一个阈值进行阈值量化处理.
这一步是最关键的一步,主要体现在阈值的选择与量化处理的过程,在每层阈值的选择上matlab提供了很多自适应的方法, 这里不一一介绍,量化处理方法主要有硬阈值量化与软阈值量化。下图是二者的区别:
上面左图是硬阈值量化,右图是软阈值量化。采用两种不同的方法,达到的效果是,硬阈值方法可以很好地保留信号边缘等局部特征,软阈值处理相对要平滑,但会造成边缘模糊等失真现象。
信号的小波重构。根据小波分解的第 N层的低频系数和经过量化处理后的第1层到第N 层的高频系数,进行信号的小波重构。
小波阀值去噪的基本问题包括三个方面:小波基的选择,阀值的选择,阀值函数的选择。
(1)小波基的选择:通常我们希望所选取的小波满足以下条件:正交性、高消失矩、紧支性、对称性或反对称性。但事实上具有上述性质的小波是不可能存在的,因为小波是对称或反对称的只有Haar小波,并且高消失矩与紧支性是一对矛盾,所以在应用的时候一般选取具有紧支的小波以及根据信号的特征来选取较为合适的小波。
(2)阀值的选择:直接影响去噪效果的一个重要因素就是阀值的选取,不同的阀值选取将有不同的去噪效果。目前主要有通用阀值(VisuShrink)、SureShrink阀值、Minimax阀值、BayesShrink阀值等。
(3)阀值函数的选择:阀值函数是修正小波系数的规则,不同的反之函数体现了不同的处理小波系数的策略。最常用的阀值函数有两种:一种是硬阀值函数,另一种是软阀值函数。还有一种介于软、硬阀值函数之间的Garrote函数。
另外,对于去噪效果好坏的评价,常用信号的信噪比(SNR)与估计信号同原始信号的均方根误差(RMSE)来判断。
参考:
小波阀值去噪法基础 http://blog.sina.com.cn/s/blog_4d7c97a00101cib3.html
在python中使用小波分析进行阈值去噪声,使用pywt.threshold函数
#coding=gbk
#使用小波分析进行阈值去噪声,使用pywt.threshold
import pywt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
data = np.linspace(1, 10, 10)
print(data)
# [ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
# pywt.threshold(data, value, mode, substitute) mode 模式有4种,soft, hard, greater, less; substitute是替换值
data_soft = pywt.threshold(data=data, value=6, mode='soft', substitute=12)
print(data_soft)
# [12. 12. 12. 12. 12. 0. 1. 2. 3. 4.] 将小于6 的值设置为12, 大于等于6 的值全部减去6
data_hard = pywt.threshold(data=data, value=6, mode='hard', substitute=12)
print(data_hard)
# [12. 12. 12. 12. 12. 6. 7. 8. 9. 10.] 将小于6 的值设置为12, 其余的值不变
data_greater = pywt.threshold(data, 6, 'greater', 12)
print(data_greater)
# [12. 12. 12. 12. 12. 6. 7. 8. 9. 10.] 将小于6 的值设置为12,大于等于阈值的值不变化
data_less = pywt.threshold(data, 6, 'less', 12)
print(data_less)
# [ 1. 2. 3. 4. 5. 6. 12. 12. 12. 12.] 将大于6 的值设置为12, 小于等于阈值的值不变
1、小波变换的理解
傅里叶变换——短时傅里叶变换——小波变换。
参考文献:以下两篇参考资料讲述得十分清楚,有助于理解小波变换。
但具体的数学角度阐述,请参考其他资料。
(1)知乎专栏:形象易懂讲解算法I——小波变换
https://zhuanlan.zhihu.com/p/22450818
(2)知乎专栏:傅里叶分析之掐死教程。
https://zhuanlan.zhihu.com/p/19763358
2、小波包分解
小波包是为了克服小波分解在高频段的频率分辨率较差,而在低频段的时间分辨率较差的问题的基础上而提出的。
它是一种更精细的信号分析的方法,提高了信号的时域分辨率。
下面是两者的对比图:
3、能量谱
基于小波包分解提取多尺度空间能量特征的原理是把不同分解尺度上的信号能量求解出来,将这些能量值按尺度顺序排列成特征向量供识别使用。
20180510补充更新:具体计算公式如下所示,本文中未使用重构后的系数进行能量值计算,直接使用小波包分解后的系数,参考文献《基于小波包能量特征的滚动轴承故障监测方法 》。
4、Matlab代码
给出两部分代码,写成两个函数。一个是小波包分解与重构,另一个是能量谱函数。
下载地址:https://download.csdn.net/download/ckzhb/10030651
代码名称:wavelet_packetdecomposition_reconstruct
function wpt= wavelet_packetdecomposition_reconstruct( x,n,wpname )
%% 对信号进行小波包分解,得到节点的小波包系数。然后对每个节点系数进行重构。
% Decompose x at depth n with wpname wavelet packets.using Shannon entropy.
%
% x-input signal,列向量。
% n-the number of decomposition layers
% wpname-a particular wavelet.type:string.
%
%Author hubery_zhang
%Date 20170714
%%
wpt=wpdec(x,n,wpname);
% Plot wavelet packet tree (binary tree)
plot(wpt)
%% wavelet packet coefficients.default:use the front 4.
cfs0=wpcoef(wpt,[n 0]);
cfs1=wpcoef(wpt,[n 1]);
cfs2=wpcoef(wpt,[n 2]);
cfs3=wpcoef(wpt,[n 3]);
figure;
subplot(5,1,1);
plot(x);
title('原始信号');
subplot(5,1,2);
plot(cfs0);
title(['结点 ',num2str(n) ' 1',' 系数'])
subplot(5,1,3);
plot(cfs1);
title(['结点 ',num2str(n) ' 2',' 系数'])
subplot(5,1,4);
plot(cfs2);
title(['结点 ',num2str(n) ' 3',' 系数'])
subplot(5,1,5);
plot(cfs3);
title(['结点 ',num2str(n) ' 4',' 系数'])
%% reconstruct wavelet packet coefficients.
rex0=wprcoef(wpt,[n 0]);
rex1=wprcoef(wpt,[n 1]);
rex2=wprcoef(wpt,[n 2]);
rex3=wprcoef(wpt,[n 3]);
figure;
subplot(5,1,1);
plot(x);
title('原始信号');
subplot(5,1,2);
plot(rex0);
title(['重构结点 ',num2str(n) ' 1',' 系数'])
subplot(5,1,3);
plot(rex1);
title(['重构结点 ',num2str(n) ' 2',' 系数'])
subplot(5,1,4);
plot(rex2);
title(['重构结点 ',num2str(n) ' 3',' 系数'])
subplot(5,1,5);
plot(rex3);
title(['重构结点 ',num2str(n) ' 4',' 系数'])
end
代码名称:wavelet_energy_spectrum
function E = wavelet_energy_spectrum( wpt,n )
%% 计算每一层每一个节点的能量
% wpt-wavelet packet tree
% n-第n层能量
%
% Author hubery_zhang
% Date 20170714
%%
% 求第n层第i个节点的系数
E(1:2^n )=0;
for i=1:2^n
E(i) = norm(wpcoef(wpt,[n,i-1]),2)^2; %20180604更新 原代码:E(i) = norm(wpcoef(wpt,[n,i-1]),2)
end
%求每个节点的概率
E_total=sum(E);
for i=1:2^n
p_node(i)= 100*E(i)/E_total;
end
% E = wenergy(wpt); only get the last layer
figure;
x=1:2^n;
bar(x,p_node);
title(['第',num2str(n),'层']);
axis([0 2^n 0 100]);
xlabel('结点');
ylabel('能量百分比/%');
for j=1:2^n
text(x(j),p_node(i),num2str(p_node(j),'%0.2f'),...
'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
end
end
一维、二维离散卷积的计算方法:
一、定义
离散信号f(n),g(n)的定义如下:
N-----为信号f(n)的长度
s(n)----为卷积结果序列,长度为len(f(n))+len(g(n))-1
以3个元素的信号为例:
f(n) = [1 2 3]; g(n) = [2 3 1];
s(0) = f(0)g(0-0) + f(1)g(0-1)+f(2)g(0-2) = 1*2 + 2*0 + 3*0 =2
s(1) = f(0)g(1-0) + f(1)g(1-1) + f(2)g(1-2) = 1*3 + 2*2 + 3*0 = 7
s(2) = f(0)g(2-0) + f(1)g(2-1) + f(2)g(2-2) =1*1 + 2*3 + 3*2=13
s(3) = f(0)g(3-0) + f(1)g(3-1) + f(2)g(3-2) =1*0 + 2*1 + 3*3=11
s(4) = f(0)g(4-0) + f(1)g(4-1) + f(2)g(4-2) =1*0 + 2*0 + 3*1=3
最终结果为:
s(n) = [2 7 13 11 3]
上述计算图示如下:
在数学里我们知道f(-x)的图像是f(x)对y轴的反转
g(-m)就是把g(m)的序列反转,g(n-m)的意义是把g(-m)平移的n点:
如上图g(m)在信号处理中通常叫做滤波器或掩码,卷积相当于掩码g(m)反转后在信号f(n)上平移求和。Matlab计算卷积的函数为conv,
>> A = [1 2 3];
B = [2,3,1];
convD = conv(A,B)
convD =
2 7 13 11 3
相应的二维卷积定义如下:
python小波去噪的方法_小波去噪基本概念相关推荐
- 小程序模板网站平台_小程序模板平台哪个好
小程序模板网站平台_小程序模板平台哪个好?分享一个微信小程序模板平台,超60个行业的小程序模板免费使用,页面内容丰富样式多样的,小程序界面模板. 微信小程序模板网站平台 微信小程序模板平台的存在,就是 ...
- python小波分解与重构_小波分解和重构
小波变换能够很好地表征一大类以低频信息为主要成分的信号, 小波包变换可以对高频部分提供更精细的分解 详见(http://www.cnblogs.com/welen/articles/5667217.h ...
- 单片机c语言小波阈值降噪,小波阈值去噪的基本原理_小波去噪阈值如何选取
小波阈值去噪的基本原理 小波阈值去噪的基本思想是先设置一个临界阈值λ,若小波系数小于λ,认为该系数主要由噪声引起,去除这部分系数;若小波系数大于λ,则认为此系数主要是由信号引起,保留这部分系数,然后对 ...
- python少儿编程讲师笔试题_小码王教育儿童编程教师面试:做笔试题(填空题和编程题,填空题 - 职朋职业圈...
为了帮助职业圈网友能够及时了解小码王教育的面试流程以及面试过程所涉及的面试问题,职业圈小编把刚获得的小码王教育面试经验马上编辑好,快速提供给大家,以便能够尽快帮助到有需要的人.这次面试总共花了1天.面 ...
- python分析数据差异的方法_用Python的两种方法进行方差分析
在进行数据分析时,我们往往会遇到要对某个变量的影响因素进行分析的情况,而影响一事物的因素往往是很多的.比如在化工生产中,有温度.压力.剂量.反应时间等因素.每一因素的改变都有可能影响产品的数量和质量. ...
- python的pygame库使用方法_[宜配屋]听图阁
使用python pygame库实现一个双人弹球小游戏,两人分别控制一个左右移动的挡板用来拦截小球,小球会在两板间不停弹跳,拦截失败的一方输掉游戏,规则类似于简化版的乒乓球. 因为是第一次用pygam ...
- 公众号跳转小程序首次没有数据_小程序如何从“0”开始运营,变成获客神器...
随着发展,越来越多的企业都感觉到获客难,获客成本高.而小程序的诞生恰恰解决了这些问题.合理的利用小程序的功能,可以帮助商家低成本高效获客,今天我们就来谈谈具体怎么用小程序来拉新引流. 合理的利用小程序 ...
- 小程序接入h5页面_小程序与H5如何互相跳转
由于小程序官方没有提供外部H5网页直接跳转到小程序的api,所以目前只支持小程序内嵌H5,并且只有内嵌的H5才能跳回小程序 小程序跳转H5 需要用到小程序的web-view,官方文档链接 web-vi ...
- 小程序webview不全屏_小程序不在小(深度)
原标题:小程序不在小(深度) 你问:"微信小程序适合哪些行业?",回答是:"所有行业!" 你可以想一下那些做过APP的公司,不管是任何行业的公司都可以拥有属于自 ...
- python读csv最快方法_使用Python读写csv文件的三种方法
Python读写csv文件 觉得有用的话,欢迎一起讨论相互学习~Follow Me 前言 逗号分隔值(Comma-Separated Values,CSV,有时也称为字符分隔值,因为分隔字符也可以不是 ...
最新文章
- 基于架构的上网行为管理产品界面对比
- 用JAVA操作ClearCase
- LeetCode 705. Design HashSet (设计哈希集合)
- Serverless的理解
- CompareAndSwap原子操作原理
- 使用 jQuery Mobile 与 HTML5 开发 Web App (十六) —— HTML5 Web Storage
- java打jar包的方式,jar命令,maven
- mysql status lock_MySQL性能突发事件问题排查技巧
- GBDT算法原理以及实例理解(含Python代码简单实现版)
- 《GNU Emacs Lisp编程入门》读书笔记
- C# 学习笔记-Chart控件使用
- Android Studio 4.1中的模板插件
- Travis CI(持续集成)
- Golden Software BLN文件格式
- dhcp服务器响应慢,请教下大家DHCP获取慢的原因和解决思路
- 百度地图获取行政区边界坐标
- STM32H750VB程序无法下载的问题
- 巅峰战舰服务器维护中,维护公告~
- 高性价比才是王道 三大流行趋势机巅峰对垒
- CVTE2016校园大使火热招募!