雷达原理---时频分析--5.小波变换在信号降噪中的应用

信号降噪

小波分析的重要应用之一就是用于信号降噪。

一、 小波变换实现信号降噪的基本原理

含噪的一维信号模型可以表示如下：
$$s(k)=f(k)+\epsilon \cdot e(k),k=0,1,\cdots ,n-1$$
式中$s(k)$为含噪信号，$f(k)$为有用信号，$e(k)$为噪声信号。假设$e(k)$是一个高斯白噪声，通常表现为高频信号，而工程实际中$f(k)$通常为低频信号或者是一些比较平稳的信号。

因此可按如下方法进行降噪处理：首先对信号进行小波分解，由于噪声信号多包含在具有较高频率的细节中，从而可以利用门限、阈值等形式对分解所得的小波系数进行处理，然后对信号进行小波重构即可达到对信号降噪的目的。

二、 信号降噪的过程

对信号降噪实质上是抑制信号中的无用部分，增强信号中有用部分的过程。一般地，一维信号降噪的过程可分为如下3个步骤：

1. 一维信号的小波分解。选择一个小波并确定分解的层次，然后进行分解计算。
2. 小波分解高频系数的阈值量化。对各个分解尺度下的高频系数选择一个阈值进行软阈值量化处理。
3. 一维小波重构。根据小波分解的最底层低频系数和各层高频系数进行一维小波重构。

这3个步骤中，最关键的是如何选择阈值以及进行阈值量化。在某种程度上，它关系到信号降噪的质量。

三、具体分析

1. 噪声在小波分解下的特性

总体上，对于一维离散信号而言，其高频部分所影响的是小波分解的第一层细节，其低频部分影响的是小波分解的最深层和低频层。如果对一个仅由白噪声所组成的信号进行分析，则可得出结论：高频系数的幅值随着分解层次的增加而迅速地衰减，且其方差也有同样的变化趋势。

用$C(j,k)$表示噪声经过小波分解后的系数，其中$j$表示尺度，$k$表示时间。下面将噪声看作普通信号，分析它的相关性、频谱和频率分布这3个主要特征。

1. 如果分析的信号s是一个平稳、零均值的白噪声，那么它的小波分解系数是相互独立的。
2. 如果信号s是一个高斯型噪声，那么其小波分解系数是互不相关的，且服从高斯分布。
3. 如果信号s是一个平稳、有色、零均值的高斯型噪声序列，那么它的小波分解系数也是高斯序列，并且对每一个分解尺度$j$，其相应的系数也是一个平稳、有色的序列。
4. 如果信号s是一个固定的、零均值的ARMA模型，那么对每一个小波分解尺度$j$，$C(j,k)，k\in Z$也是固定的、零均值ARMA模型，且其特性取决于尺度$j$。
5. 如果s是一个一般的噪声：

• 若它的相关函数已知，则可计算系数序列$C(j,k)$和$C(j,k^{\prime })$;
• 若它的相干函数谱已知，则可计算$C(j,k)，k\in Z$的谱及尺度$j$和$j^{\prime }$的交叉谱。

2. 应用一维小波分析进行信号消噪处理

小波分析进行降噪处理一般有下述3种方法：

1. 默认阈值降噪处理。该方法利用函数ddencmp生成信号的默认阈值，然后利用函数wdencmp进行降噪处理。
2. 给定阈值降噪处理。在实际的降噪处理过程中，阈值往往可通过经验公式获得，且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可利用函数wthresh。
3. 强制降噪处理。该方法是将小波分解结构中的高频系数全部置为0，即滤掉所有高频部分，然后对信号进行小波重构。这种方法比较简单，且去噪后的信号比较平滑，但是容易丢失信号中的有用成分。

3. 消噪阈值选取规则

wthresh：对输入信号进行软阈值或硬阈值处理；

 调用格式：
  Y = wthresh(X,sorh,T);

 使用说明：

• sorh='s'时，进行软阈值处理，即把信号的绝对值与指定的阈值进行比较，小于或等于阈值的点变为0；大于阈值的点变为该点值与阈值的差。
• sorh='h'时，进行硬阈值处理，即把信号的绝对值与指定的阈值进行比较，小于或等于阈值的点变为0；大于阈值的点保持不变。
• T为阈值大小，一般来说，硬阈值比软阈值处理后的信号更加粗糙。

根据基本的噪声模型，阈值的选取T可基于函数T = thselect(X,tptr);

1. 输入参数tptr='rigrsure'，是一种基于Stein无偏似然估计（SURE）原理的自适应阈值选择。给定一个阈值t，得到它的似然估计，再将非似然最小化，就可得到所选的阈值。这是一种软件阈值估计器。
2. 输入参数tptr='sqtwolog'，是一种固定的阈值形式，它所产生的阈值为sqrt(2*log(length(X)))。
3. 输入参数tptr='heursure'，是前两种阈值的综合，所选择的是最优预测变量阈值。如果信噪比很小，而SURE估计有很大的噪声，此时就需采用这种固定的阈值形式。
4. 输入参数tptr='minimaxi'，也是一种固定的阈值选择形式，它所产生的是一个最小均方差的极值，而不是无误差。

对噪声进行小波分解时，也会产生高频系数，故一个信号的高频系数向量是有用信号和噪声信号的高频系数的叠加。由于SURE和minimaxi阈值选取规则较为保守（仅将部分系数置为零），因此在信号的高频信息有很少一部分在噪声范围内时，这两种阈值非常有用，可以将弱小的信号提取出来。其他两种阈值选取规则，在去除噪声时更为有效，但是也有可能将有用信号的高频部分当作噪声信号去除掉。

在实际应用中，大多数信号可能包含着许多尖峰或突变，而且噪声信号也并不是平稳的白噪声。对这种信号进行去噪处理时，传统的傅里叶变换完全是在频域中对信号进行分析，它不能给出信号在某个时间点上的变化情况，因此分辨不出信号在时间轴上的任何一个突变。但是小波分析能同时在时频域中对信号进行分析，所以它能有效地区分信号中的突变部分和噪声，从而实现对非平稳信号的去噪。

四、示例

1. 利用3种降噪处理方法对污染信号进行去噪处理

clear all;
close all;
clc;
%******************* 利用小波分析对污染信号进行去噪处理以恢复原始信号 *********************%
% 原始信号
load leleccum;      % 装载采集的信号 leleccum.mat
indx = 2000:3450; 
s = leleccum(indx); % 将信号中第2000到第3450个采样点赋给s
subplot(2,2,1);  plot(s); title('原始信号'); xlabel('样本序号 n'); ylabel('幅值 A');
% 用db1小波对原始信号进行3层分解并提取系数
[c,l] = wavedec(s,3,'db1');
a3 = appcoef (c,l,'db1',3);
d3 = detcoef (c,l,3);
d2 = detcoef (c,l,2);
d1 = detcoef (c,l,1);
% 对信号进行强制去噪处理
dd3 = zeros(1,length(d3));
dd2 = zeros(1,length(d2));
dd1 = zeros(1,length(d1));
c1 = [a3 dd3 dd2 dd1];
s1 = waverec(c1,l,'db1');
subplot(2,2,2);  plot(s1); title('强制去噪后的信号'); xlabel('样本序号 n'); ylabel('幅值 A');
% 用默认阈值对信号进行去噪处理
% 用ddencmp 函数获取信号的默认阈值
[thr,sorh,keepapp] = ddencmp('den','wv',s);
s2 = wdencmp ('gbl',c,l,'db1',3,thr,sorh,keepapp);
subplot(2,2,3);  plot(s2); title('默认阈值去噪后的信号'); xlabel('样本序号 n'); ylabel('幅值 A');
% 用给定的软阈值进行去噪处理
softd1 = wthresh (d1,'s',1.465);
softd2 = wthresh (d2,'s',1.823);
softd3 = wthresh (d3,'s',2.768);
c2 = [a3 softd3 softd2 softd1];
s3 = waverec(c2,l,'db1');
subplot(2,2,4);  plot(s3); title('给定软阈值去噪后的信号'); xlabel('样本序号 n'); ylabel('幅值 A');

2. 使用wden函数对信号进行自动降噪

clear all;
close all;
clc;
%******************* 利用小波分析对一个含噪的矩形波信号进行去噪处理 *********************%
% 设置信噪比和随机种子值
snr = 4;
init = 2055615866; 
% 产生原始信号sref和被高斯白噪声污染的信号s
[sref , s] = wnoise (1,11,snr,init);
% 用sym8小波对信号s进行3层分解并对信号进行降噪
xd = wden (s,'heursure','s','one',3,'sym8'); %选用SURE阈值模式
% 显示信号波形
subplot(3,1,1);  plot(sref); title('参考信号'); xlabel('样本序号 n'); ylabel('幅值 A');
subplot(3,1,2);  plot(s);    title('含噪信号'); xlabel('样本序号 n'); ylabel('幅值 A');
subplot(3,1,3);  plot(xd);   title('去噪信号'); xlabel('样本序号 n'); ylabel('幅值 A');