信号去噪算法Python 信号降噪算法

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

文章目录

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

• 一、应用背景
• 二、信号去噪

• 1、低通滤波去噪
• 2、小波分解去噪

• 2.1 Mallat金字塔算法
• 2.2 小波基的选取

• 2.2.1 小波基的种类
• 2.2.2 小波基的选取原则

• 2.3 小波分解去噪仿真

• 3、奇异值分解去噪

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

• 3.1.1 信号的矩阵重构
• 3.1.2 信号奇异值分解去噪流程

• 3.2 奇异值分解去噪仿真

• 二、信号回归

• 1、最小二乘回归

• 1.1 最小二乘回归建模
• 1.2 最小二乘回归仿真

• 2、岭回归

• 2.1 岭回归建模
• 2.2 岭回归仿真

• 3、高斯过程回归

• 3.1 高斯过程回归建模
• 3.2 高斯过程回归仿真

• 三、蒙特卡洛实验
• MATLAB代码

一、应用背景

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

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

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

原始信号	加噪后信号（$σ^{2}$
$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$

$\begin{bmatrix} y(1) & y(2) & \cdots & y(N-L+1) \\ y(2) & y(3) & \cdots & y(N-L+2) \\ \vdots & \vdots & \ddots & \vdots \\ y(L) & y(L+1) & \cdots & y(N) \\ \end{bmatrix}$

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

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

构建重构矩阵

构建重构矩阵

待去噪信号 y

重构矩阵 A

对重构矩阵 A 进行奇异值分解

以 2n 作为有效秩的阶次, 即选取 2n 个主奇异值

确定主频数量 n

对待去噪信号 y 进行傅里叶变换等操作

将其他奇异值置 0, 得到新的重构矩阵 B

将 B 中对应的元素相加后取平均得到降噪后的信号 y0

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 = \begin{bmatrix} x^6 & x^5 & x^4 & x^3 & x^2 & x & 1 \end{bmatrix}$

系数 $\omega$ 为：
$\omega = \begin{bmatrix} \omega_6 \\ \omega_5 \\ \omega_4 \\ \omega_3 \\ \omega_2 \\ \omega_1 \\ \omega_0 \end{bmatrix}$

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

1.2 最小二乘回归仿真

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

2、岭回归

2.1 岭回归建模

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

$\mathop {\min }\limits_\omega E\left[ {{{\left( {Y - \sum\limits_{k = 1}^N {{\omega_k} {X_k}} } \right)}^2}} \right]$

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

$\mathop {\min }\limits_\omega E \left[ {{{\left( {Y - \sum\limits_{k = 1}^N {{\omega_k}{X_k}} } \right)}^2} + \lambda \sum\limits_{k = 1}^N {\omega_k^2} } \right]$

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

$\begin{aligned} L\left( {\omega ,\lambda } \right) &= \sum\limits_{k = 1}^N {{{\left( {{Y_k} - {\omega_k}{X_k}} \right)}^2} + \lambda \cdot \left\| \omega \right\|_2^2} \\ &= {\left( {Y - X\omega } \right)^T}\left( {Y - X\omega } \right){\text{ + }}\lambda \cdot \left\| \omega \right\|_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{aligned}$

求导有：

${\nabla_\omega }L\left( {\omega ,\lambda } \right) = - {X^T}Y - {X^T}Y + 2{X^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{aligned} E\left( {{Y_{n + 1}}\left| {{Y_1},...,{Y_n}} \right.} \right) &= {\Sigma _{Y,{Y_{n + 1}}}} \cdot \Sigma _{YY}^{ - 1} \cdot E\left( Y \right) \\ Var\left( {{Y_{n + 1}}\left| {{Y_1},...,{Y_n}} \right.} \right) &= {\Sigma _{{Y_{n + 1}}}} - {\Sigma _{Y,{Y_{n + 1}}}} \cdot \Sigma _Y^{ - 1} \cdot {\Sigma _{{Y_{n + 1}},Y}} \end{aligned}$
所以可得高斯过程回归（GPR）
$\begin{bmatrix} Y \\ Y_{n+1} \end{bmatrix} \sim N\left( {0, \begin{bmatrix} \Sigma_Y & \Sigma_{Y,Y_{n+1}} \\ \Sigma_{Y_{n+1},Y} & \Sigma_{Y_{n+1}} \end{bmatrix} } \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('高斯过程回归')