常用信号去噪与信号回归方法的原理及MATLAB实现

常用信号去噪与回归方法的原理及MATLAB实现

一、应用背景

  “信号去噪”与“信号回归”是信号处理的基本技术。本博客针对一类特殊的问题（峰值缓变函数优化），对这两类技术展开讨论。
  峰值缓变函数是指最值附近的函数值变化十分平稳，当存在噪声干扰时，会导致优化算法的性能十分不稳定。因此，对于这类函数（信号）为了充分发挥优化算法的性能，可在前期进行“去噪”或“回归”处理。

例如：
$$f\left(x\right)=\frac{1}{1+{\left(x-1\right)}^6\text{ + }0.02x}$$

该函数（信号）及其加噪后的最优解如下所示：

原始信号	加噪后信号（${\sigma }^2$ =0.02）
$x_{opt}$=0.68	$x_{opt}$=1.05

接下来，便针对如何从含噪信号最大程度恢复最优解展开讨论与实验。

---

二、信号去噪

1、低通滤波去噪

  由于目标函数在峰值附近的缓变的，故峰值附近的频谱以低频分量为主；而噪声往往分布在高频区域。故直接进行低通滤波可在一定程度上实现去噪，而几乎不会改变峰值的特性。

	原始信号	含噪信号	去噪信号
信号
频谱
最优解	$x_{opt}$=0.68	$x_{opt}$=1.05	$x_{opt}$=0.66

2、小波分解去噪

2.1 Mallat金字塔算法

  Mallat金字塔算法，利用“双尺度关系”，为小波变换赋予了“多分辨率表示”的物理意义。内层分解表示信号的基本（大尺度）信息，而越往外层，分解提供越细节（小尺度）的信息。

2.2 小波基的选取

2.2.1 小波基的种类

MATLAB小波基函数（MathWorks）
MATLAB中的小波基介绍（CSDN）

2.2.2 小波基的选取原则

① 正交性
② 支撑区间
  支撑区间长度是指随 时间/频率 趋于 ∞ 时，小波函数收敛到零所需的长度。
  “紧支撑”是指对于函数 f(x)，如果自变量 x 在 0 附近的取值范围内，f(x) 能取到值；而在此之外，f(x) 取值为 0，那么这个函数 f(x) 就是紧支撑函数，而这个 0 附近的取值范围就叫做紧支撑集。总结为一句话就是“除在一个很小的区域外，函数取值为零，即函数有速降性”。
③ 对称性
  对称性使小波具有线性相移，这对于图像处理具有重要意义。
④ 正则性
  正则性描述小波基函数的平滑程度

2.3 小波分解去噪仿真

这里选取 7 层小波分解，内3层利用软阈值（0.014）处理，外4层利用硬阈值置0；

	原始信号	含噪信号	去噪信号
信号
最优解	$x_{opt}=0.68$	$x_{opt}=1.05$	$x_{opt}=0.71$

3、奇异值分解去噪

3.1 信号奇异值分解去噪原理

  前面所述的傅里叶变换与小波变换，都是从信号能量入手，通过对信号进行某种分解，实现去噪。
  而借助奇异值分解（SVD）我们可以从信号结构层面进行去噪。然而，奇异值分解是针对矩阵的，这里待去噪信号仅有一个一维信号，故我们首先需要构造信号的重构矩阵。

3.1.1 信号的矩阵重构

  对于信号 $y(i),i=1,2,...,N$，基于相空间重构理论，可以由其构造重构矩阵$A$

$$\left[\begin{array}{cccc}y(1)& y(2)& \cdots & y(N-L+1)\\ y(2)& y(3)& \cdots & y(N-L+2)\\ ⋮& ⋮& \ddots & ⋮\\ y(L)& y(L+1)& \cdots & y(N)\end{array}\right]$$

其中，不同的 L 选取，对应着不同尺寸的重构矩阵，一般 L 取为信号长度的一半。

3.1.2 信号奇异值分解去噪流程

3.2 奇异值分解去噪仿真

	原始信号	含噪信号	去噪信号
信号
最优解	$x_{opt}=0.68$	$x_{opt}=1.05$	$x_{opt}=0.71$

---

二、信号回归

  信号去噪类方法通过直接去除噪声达到恢复信号的目的。与此对应的，信号回归类方法，通过使恢复信号达到最小均方误差等约束，直接对信号进行重建。

1、最小二乘回归

  最小二乘法 (LS) 能提供对信号的最小方差无偏估计，在线性模型中，最小二乘方法是最优的，其能够达到CRLB。

1.1 最小二乘回归建模

  最小二乘针对的是线性模型，而这里的待重建信号为非线性信号。因此，我们只能通过对信号进行若干阶次的多项式拟合来重构信号。
  这里，我们采用 6 次多项式进行（因为我们有先验知识信号是 6 阶），可建立如下模型：
$$y=X\omega$$
其中，y 为含噪信号，X为：

$$X=\left[\begin{array}{ccccccc}{x}^6& x^5& x^4& x^3& x^2& x& 1\end{array}\right]$$

系数 $\omega$ 为：
$$\omega =\left[\begin{array}{c}{\omega }_6\\ {\omega }_5\\ {\omega }_4\\ {\omega }_3\\ {\omega }_2\\ {\omega }_1\\ {\omega }_0\end{array}\right]$$

容易得到，最小二乘解为：
$${\omega }_{LS}=X^{+}y=(X^{H}X{)}^{-1}{X}^{H}y$$

1.2 最小二乘回归仿真

	原始信号	含噪信号	回归信号
信号
最优解	$x_{opt}=0.68$	$x_{opt}=1.05$	$x_{opt}=0.65$

2、岭回归

2.1 岭回归建模

  虽然最小二乘是线性模型的最优解，但是其非常容易受到奇异数据或病态数据的干扰。特别是当信号噪声较大时，回归结果极易受到噪声干扰。
  这里，为了使得最小二乘适应强噪声情况，我们对最小二乘模型

$$\underset{\omega }{\min }E\left[{\left(Y-\sum _{k=1}^N{\omega }_k{X}_k\right)}^2\right]$$

施加L2正则化（岭回归）：

$$\underset{\omega }{\min }E\left[{\left(Y-\sum _{k=1}^N{\omega }_k{X}_k\right)}^2+\lambda \sum _{k=1}^N{\omega }_k^2\right]$$

下面求解岭回归的最优解，由于该优化为凸优化，故直接由一阶必要条件，求导即可得到全局最优解。

$$\begin{array}{rl}L\left(\omega ,\lambda \right)& =\sum _{k=1}^N{\left(Y_k-{\omega }_k{X}_k\right)}^2+\lambda \cdot {\Vert \omega \Vert }_2^2\\ & ={\left(Y-X\omega \right)}^T\left(Y-X\omega \right)\text{ + }\lambda \cdot {\Vert \omega \Vert }_2^2\\ & =Y^{T}Y-{\omega }^T{X}^{T}Y-Y^{T}X\omega +{\omega }^T{X}^{T}X\omega +\lambda {\omega }^T\omega \end{array}$$

求导有：

$${\nabla }_{\omega }L\left(\omega ,\lambda \right)=-X^{T}Y-X^{T}Y+2X^{T}X\omega +2\lambda \omega$$

令

$${\nabla }_{\omega }L\left(\omega ,\lambda \right)=0$$

可得

$${\omega }_{Ridge}={\left(X^{T}X+\lambda I\right)}^{-1}{X}^{T}Y$$

可见，只需调节$\lambda$，便可避免对于噪声的“过拟合”。

2.2 岭回归仿真

	原始信号	含噪信号	回归信号
信号
最优解	$x_{opt}=0.68$	$x_{opt}=1.05$	$x_{opt}=0.68$

3、高斯过程回归

3.1 高斯过程回归建模

  若$Y_1,Y_2,...,Y_n$服从高斯分布，即它们都是高斯过程 $Y(t)$ 的采样，则由高斯过程的条件期望，对于 $Y=(Y_1,Y_2,...,Y_n)$ ，预测$Y_{n+1}$有：
$$\begin{array}{rl}E\left(Y_{n+1}|Y_1,...,Y_n\right)& ={\Sigma }_{Y,Y_{n+1}}\cdot {\Sigma }_{YY}^{-1}\cdot E\left(Y\right)\\ Var\left(Y_{n+1}|Y_1,...,Y_n\right)& ={\Sigma }_{Y_{n+1}}-{\Sigma }_{Y,Y_{n+1}}\cdot {\Sigma }_Y^{-1}\cdot {\Sigma }_{Y_{n+1},Y}\end{array}$$
所以可得高斯过程回归（GPR）
$$\left[\begin{array}{c}Y\\ Y_{n+1}\end{array}\right]\sim N\left(0,\left[\begin{array}{cc}{\Sigma }_Y& {\Sigma }_{Y,Y_{n+1}}\\ {\Sigma }_{Y_{n+1},Y}& {\Sigma }_{Y_{n+1}}\end{array}\right]\right)$$
从而可利用历史数据训练得到相关矩阵，再由相关矩阵预测$Y_{n+1}$的均值与方差。逐步迭代，完成对完整信号的回归。

3.2 高斯过程回归仿真

  尽管这里的信号 y(t) 并不是高斯过程，但其噪声是高斯过程。另外，值得一提的是高斯过程回归总是能取得“意想不到”的结果。

	原始信号	含噪信号	回归信号
信号
最优解	$x_{opt}=0.68$	$x_{opt}=1.05$	$x_{opt}=0.70$

---

三、蒙特卡洛实验

对于前面几节所介绍的方法，取噪声方差 0.02，进行蒙特卡洛实验，得到它们对于信号恢复能力的对比，见下表。

	平均MSE	平均最优解	最优解平均偏差
原始信号	0	69	0
含噪信号	0.0198	80.46	14.92
低通滤波去噪	0.0098	67.36	5.42
小波分解去噪	0.0059	72.73	6.25
奇异值分解去噪	0.0018	70.47	3.41
最小二乘回归	0.0043	66.62	2.38
岭回归	0.0033	69.04	0.42
高斯过程回归	0.0011	70.31	7.29

---

MATLAB代码

上述全部的内容的MATLAB代码如下。
（最后，如有任何问题，欢迎讨论与指导。）

% 峰值缓变函数最优化（去噪方法、回归方法）
% CopyRight @ TSENG ChihYuan
%
clear,clc,close all
%% 构造峰值缓变函数
x = 0 : 0.01 : 2;
%y = 1 ./ ( 1 + ( (x-1).^6) ); % 巴特沃斯低通滤波函数（峰值缓变）
y = 1 ./ ( 1 + ( (x-1).^6 + 0.02 * x ) );
% 最值求解
pos = find( y == max(y) );%寻找最大值
disp([ '原始曲线 : ','最优解=' , num2str( x(pos) ) ,'， 最优函数值=' , num2str( y(pos) ) ])
%% 加噪（高斯白噪声）
sigma = 0.01; %噪声方差
noise = normrnd(0,sigma,1,length(x));
yn = y + noise ;
% 最值求解
pos2 = find( yn == max(yn) );%寻找最大值
disp([ '加噪曲线 : ','最优解=' , num2str( x(pos2) ) ,'， 最优函数值=' , num2str( y(pos2) ) ])
%% %%%%%%%%%%%%%%%%%%%%%%%%% 峰值缓变函数最优化（去噪思路）%%%%%%%%%%%%%%%%%%%%%%%%
disp(['------- 去噪方法 -------'])
%% 低通滤波去噪
f2 = fft(yn);
f2(7:193) = 0; %这里选取的主频个数为7（后面的方法和这里保持一致）
y3 = ifft(f2);
% 最值求解
pos3 = find( y3 == max(y3) );%寻找最大值
disp([ '低通滤波去噪 : ','最优解=' , num2str( x(pos3) ) ,'， 最优函数值=' , num2str( y(pos3) ) ])
%% 小波分解去噪
% https://ww2.mathworks.cn/help/wavelet/ref/wfilters.html
[c,l] = wavedec(yn,7,'coif5'); %Mallat小波分解（分解层数为7层，小波基为coif5）
ca11=appcoef(c,l,'coif5',7); %获取低频信号
%获取高频细节
cd1=detcoef(c,l,1);
cd2=detcoef(c,l,2); 
cd3=detcoef(c,l,3);
cd4=detcoef(c,l,4);
cd5=detcoef(c,l,5);
cd6=detcoef(c,l,6);
cd7=detcoef(c,l,7);
% 1-3层置0,4-7层用软阈值函数处理
sd1=zeros(1,length(cd1));
sd2=zeros(1,length(cd2)); 
sd3=zeros(1,length(cd3));
sd4=wthresh(cd4,'s',0.014);
sd5=wthresh(cd5,'s',0.014);
sd6=wthresh(cd6,'s',0.014);
sd7=wthresh(cd7,'s',0.014);
c2=[ca11,sd7,sd6,sd5,sd4,sd3,sd2,sd1];
y4=waverec(c2,l,'coif5'); %小波重构 
% 最值求解
pos4 = find( y4 == max(y4) );%寻找最大值
disp([ '小波分解去噪 : ','最优解=' , num2str( x(pos4) ) ,'， 最优函数值=' , num2str( y(pos4) ) ])
%% 奇异值分解去噪
N = length(yn);%信号长度
A = CorMat(yn,floor(N/2));%由信号构造相关矩阵
[U,S,V] = svd(A); %SVD
S( 5:floor(N/2),5:floor(N/2) ) = 0; %这里主元数选5
B = U * S * V';
y5 = CorMat2(B);
% 最值求解
pos5 = find( y5 == max(y5) );%寻找最大值
disp([ '奇异值分解去噪 : ','最优解=' , num2str( x(pos5) ) ,'， 最优函数值=' , num2str( y(pos5) ) ])


%% %%%%%%%%%%%%%%%%%%%%%%%%% 峰值缓变函数最优化（回归思路）%%%%%%%%%%%%%%%%%%%%%%%%
disp(['------- 回归方法 -------'])
%% 最小二乘
% y = X * c 其中，yn为 n×1 向量，X为 n×7矩阵 ，c为 7×1向量
x = x'; yn = yn'; 
X = [ x.^6 , x.^5 , x.^4 , x.^3 , x.^2 , x.^1 , x.^0 ]; %这里采用6多项式（因为原函数是6次）
c = pinv(X) * yn; %LS
yLS = X * c;
% 最值求解
pos13 = find( yLS == max(yLS) );%寻找最大值
disp([ '最小二乘回归 : ','最优解=' , num2str( x(pos13) ) ,'， 最优函数值=' , num2str( y(pos13) ) ])
%% 岭回归
lamda = 0.001;
X = [ x.^6 , x.^5 , x.^4 , x.^3 , x.^2 , x.^1 , x.^0 ]; %这里采用6多项式（因为原函数是6次）
E = eye(size(X'*X));
w = (X'*X + lamda * E) \ ( X'*yn ); %岭回归
yRidge = X * w;
% 最值求解
pos14 = find( yRidge == max(yRidge) );%寻找最大值
disp([ '岭回归 : ','最优解=' , num2str( x(pos14) ) ,'， 最优函数值=' , num2str( y(pos14) ) ])
%% 高斯过程回归
% fitrgp()的说明页面 - https://ww2.mathworks.cn/help/stats/fitrgp.html#butnn96-2
gprMdl = fitrgp(x,yn);
%gprMdl = fitrgp(x,yn,'Basis','pureQuadratic','FitMethod','exact','PredictMethod','exact');
yGauss = resubPredict(gprMdl);
% 最值求解
pos15 = find( yGauss == max(yGauss) );%寻找最大值
disp([ '高斯过程回归 : ','最优解=' , num2str( x(pos15) ) ,'， 最优函数值=' , num2str( y(pos15) ) ])

%% 绘图
figure,plot(x,y,'.-'),title('原始峰值缓变函数')
figure,plot(x,yn,'.-'),title('原函数加噪')
figure,plot(x,abs(y3),'.-'),title('低通滤波去噪')
figure,plot(x,y4,'.-'),title('小波分解去噪')
figure,plot(x,y5,'.-'),title('奇异值分解去噪')
figure,plot(x,yLS,'.-'),title('最小二乘回归')
figure,plot(x,yRidge,'.-'),title('岭回归')
figure,plot(x,yGauss,'.-'),title('高斯过程回归')