小波(Wavelet)这一术语,顾名思义,“小波”就是小的波形。所谓“小”是指它具有衰减性;而称之为“波”则是指它的波动性,其振幅正负相间的震荡形式。与Fourier变换相比,小波变换是时间(空间)频率的局部化分析,它通过伸缩平移运算对信号(函数)逐步进行多尺度细化,最终达到高频处时间细分,低频处频率细分,能自动适应时频信号分析的要求,从而可聚焦到信号的任意细节,解决了Fourier变换的困难问题,成为继Fourier变换以来在科学方法上的重大突破。有人把小波变换称为“数学显微镜”。
小波变换的概念是由法国从事石油信号处理的工程师J.Morlet在1974年首先提出的,通过物理的直观和信号处理的实际经验的需要建立了反演公式,当时未能得到数学家的认可。正如1807年法国的热学工程师J.B.J.Fourier提出任一函数都能展开成三角函数的无穷级数的创新概念未能得到著名数学家J.L.Lagrange,P.S.Laplace以及A.M.Legendre的认可一样。幸运的是,早在七十年代,A.Calderon表示定理的发现、Hardy空间的原子分解和无条件基的深入研究为小波变换的诞生做了理论上的准备,而且J.O.Stromberg还构造了历史上非常类似于当前的小波基;1986年著名数学家Y.Meyer偶然构造出一个真正的小波基,并与S.Mallat合作建立了构造小波基的统一方法加多尺度分析之后,小波分析才开始蓬勃发展起来,其中比利时女数学家I.Daubechies撰写的《小波十讲(Ten Lectures on Wavelets)》对小波的普及起了重要的推动作用。它与Fourier变换、窗口Fourier变换(Gabor变换)相比,这是一个时间和频率的局域变换,因而能有效的从信号中提取信息,通过伸缩和平移等运算功能对函数或信号进行多尺度细化分析(Multiscale Analysis),解决了Fourier变换不能解决的许多困难问题,从而小波变化被誉为“数学显微镜”,它是调和分析发展史上里程碑式的进展。
1简介
所谓的小波就是指微小的具有波动性的波形,从信号学角度来看,是一个信号滤波的过程。目前,小波去噪已然是一种喜闻乐见的噪声处理方式。
本文使用matlab进行小波降噪,采用软阈值、硬阈值、和固定阈值三种方式,对一维噪声数据进行降噪处理,并输出误差SNR,RMSE,以及降噪前后的噪声图像 。
2 操作步骤
主要有三个步骤
对被噪声污染的信号进行小波变换处理;
通过对变换后的小波系数进行处理,去除其中的噪声
小波逆变换,得到去噪后的信号
3 代码
%%初始化程序
clear,clc
t1=clock;
%% 载入噪声信号数据,数据为.mat格式,并且和程序放置在同一个文件夹下
load('filename.mat');%matrix
YSJ= filename;
%% 数据预处理,数据可能是存储在矩阵或者是EXCEL中的二维数据,衔接为一维的,如果数据是一维数据,此步骤也不会影响数据
[c,l]=size(YSJ);
Y=[];
for i=1:c
Y=[Y,YSJ(i,:)];
end
[c1,l1]=size(Y);
X=[1:l1];
%% 绘制噪声信号图像
figure(1);
plot(X,Y);
xlabel('横坐标');
ylabel('纵坐标');
title('原始信号');
%% 硬阈值处理
lev=3;
xd=wden(Y,'heursure','h','one',lev,'db4');%硬阈值去噪处理后的信号序列
figure(2)
plot(X,xd)
xlabel('横坐标');
ylabel('纵坐标');
title('硬阈值去噪处理')
set(gcf,'Color',[1 1 1])
%% 软阈值处理
lev=3;
xs=wden(Y,'heursure','s','one',lev,'db4');%软阈值去噪处理后的信号序列
figure(3)
plot(X,xs)
xlabel('横坐标');
ylabel('纵坐标');
title('软阈值去噪处理')
set(gcf,'Color',[1 1 1])
%% 固定阈值后的去噪处理
lev=3;
xz=wden(Y,'sqtwolog','s','sln',lev,'db4');%固定阈值去噪处理后的信号序列
figure(4)
plot(X,xz);
xlabel('横坐标');
ylabel('纵坐标');
title('固定阈值后的去噪处理')
set(gcf,'Color',[1 1 1])
%% 计算信噪比SNR
Psig=sum(Y*Y')/l1;
Pnoi1=sum((Y-xd)*(Y-xd)')/l1;
Pnoi2=sum((Y-xs)*(Y-xs)')/l1;
Pnoi3=sum((Y-xz)*(Y-xz)')/l1;
SNR1=10*log10(Psig/Pnoi1);
SNR2=10*log10(Psig/Pnoi2);
SNR3=10*log10(Psig/Pnoi3);
%% 计算均方根误差RMSE
RMSE1=sqrt(Pnoi1);
RMSE2=sqrt(Pnoi2);
RMSE3=sqrt(Pnoi3);
%% 输出结果
disp('-------------三种阈值设定方式的降噪处理结果---------------');
disp(['硬阈值去噪处理的SNR=',num2str(SNR1),',RMSE=',num2str(RMSE1)]);
disp(['软阈值去噪处理的SNR=',num2str(SNR2),',RMSE=',num2str(RMSE2)]);
disp(['固定阈值后的去噪处理SNR=',num2str(SNR3),',RMSE=',num2str(RMSE3)]);
t2=clock;
tim=etime(t2,t1);
disp(['------------------运行耗时',num2str(tim),'秒-------------------'])
4 结果
1. 噪声信号图像
2. 硬阈值去噪图像
3. 软阈值去噪图像
4. 固定阈值去噪图像
5. 去噪结果的误差对比