常用信号去噪与回归方法的原理及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(x)=11+(x−1)6+ 0.02xf\left( x \right) = \frac{1}{{1 + {{\left( {x - 1} \right)}^6}{\text{ + }}0.02x}} f(x)=1+(x−1)6 + 0.02x1

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

原始信号	加噪后信号（σ2σ^{2}σ2 =0.02）

xoptx_{opt}xopt=0.68	xoptx_{opt}xopt=1.05

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

二、信号去噪

1、低通滤波去噪

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

	原始信号	含噪信号	去噪信号
信号
频谱
最优解	xoptx_{opt}xopt=0.68	xoptx_{opt}xopt=1.05	xoptx_{opt}xopt=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；

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

3、奇异值分解去噪

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

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

3.1.1 信号的矩阵重构

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

[y(1)y(2)⋯y(N−L+1)y(2)y(3)⋯y(N−L+2)⋮⋮⋱⋮y(L)y(L+1)⋯y(N)]\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} ⎣⎢⎢⎢⎡y(1)y(2)⋮y(L)y(2)y(3)⋮y(L+1)⋯⋯⋱⋯y(N−L+1)y(N−L+2)⋮y(N)⎦⎥⎥⎥⎤

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

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

构建重构矩阵

待去噪信号 y

重构矩阵 A

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

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

确定主频数量 n

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

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

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

3.2 奇异值分解去噪仿真

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

二、信号回归

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

1、最小二乘回归

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

1.1 最小二乘回归建模

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

X=[x6x5x4x3x2x1]X = \begin{bmatrix} x^6 & x^5 & x^4 & x^3 & x^2 & x & 1 \end{bmatrix} X=[x6x5x4x3x2x1]

系数 ω\omegaω 为：
ω=[ω6ω5ω4ω3ω2ω1ω0]\omega = \begin{bmatrix} \omega_6 \\ \omega_5 \\ \omega_4 \\ \omega_3 \\ \omega_2 \\ \omega_1 \\ \omega_0 \end{bmatrix} ω=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡ω6ω5ω4ω3ω2ω1ω0⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤

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

1.2 最小二乘回归仿真

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

2、岭回归

2.1 岭回归建模

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

min⁡ωE[(Y−∑k=1NωkXk)2]\mathop {\min }\limits_\omega E\left[ {{{\left( {Y - \sum\limits_{k = 1}^N {{\omega_k} {X_k}} } \right)}^2}} \right] ωminE⎣⎡(Y−k=1∑NωkXk)2⎦⎤

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

min⁡ωE[(Y−∑k=1NωkXk)2+λ∑k=1Nωk2]\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] ωminE⎣⎡(Y−k=1∑NωkXk)2+λk=1∑Nωk2⎦⎤

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

L(ω,λ)=∑k=1N(Yk−ωkXk)2+λ⋅∥ω∥22=(Y−Xω)T(Y−Xω)+ λ⋅∥ω∥22=YTY−ωTXTY−YTXω+ωTXTXω+λωTω\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} L(ω,λ)=k=1∑N(Yk−ωkXk)2+λ⋅∥ω∥22=(Y−Xω)T(Y−Xω) + λ⋅∥ω∥22=YTY−ωTXTY−YTXω+ωTXTXω+λωTω

求导有：

∇ωL(ω,λ)=−XTY−XTY+2XTXω+2λω{\nabla_\omega }L\left( {\omega ,\lambda } \right) = - {X^T}Y - {X^T}Y + 2{X^T}X\omega + 2\lambda \omega ∇ωL(ω,λ)=−XTY−XTY+2XTXω+2λω

令

∇ωL(ω,λ)=0{\nabla_\omega }L\left( {\omega ,\lambda } \right) = 0 ∇ωL(ω,λ)=0

可得

ωRidge=(XTX+λI)−1XTY{\omega_{Ridge}} = {\left( {{X^T}X + \lambda I} \right)^{ - 1}}{X^T}Y ωRidge=(XTX+λI)−1XTY

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

2.2 岭回归仿真

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

3、高斯过程回归

3.1 高斯过程回归建模

若Y1,Y2,...,YnY_1,Y_2,...,Y_nY1,Y2,...,Yn服从高斯分布，即它们都是高斯过程 Y(t)Y(t)Y(t) 的采样，则由高斯过程的条件期望，对于 Y=(Y1,Y2,...,Yn)Y=(Y_1,Y_2,...,Y_n)Y=(Y1,Y2,...,Yn) ，预测Yn+1Y_{n+1}Yn+1有：
E(Yn+1∣Y1,...,Yn)=ΣY,Yn+1⋅ΣYY−1⋅E(Y)Var(Yn+1∣Y1,...,Yn)=ΣYn+1−ΣY,Yn+1⋅ΣY−1⋅ΣYn+1,Y\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} E(Yn+1∣Y1,...,Yn)Var(Yn+1∣Y1,...,Yn)=ΣY,Yn+1⋅ΣYY−1⋅E(Y)=ΣYn+1−ΣY,Yn+1⋅ΣY−1⋅ΣYn+1,Y
所以可得高斯过程回归（GPR）
[YYn+1]∼N(0,[ΣYΣY,Yn+1ΣYn+1,YΣYn+1])\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) [YYn+1]∼N(0,[ΣYΣYn+1,YΣY,Yn+1ΣYn+1])
从而可利用历史数据训练得到相关矩阵，再由相关矩阵预测Yn+1Y_{n+1}Yn+1的均值与方差。逐步迭代，完成对完整信号的回归。

3.2 高斯过程回归仿真

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

	原始信号	含噪信号	回归信号
信号
最优解	xopt=0.68x_{opt}=0.68xopt=0.68	xopt=1.05x_{opt}=1.05xopt=1.05	xopt=0.70x_{opt}=0.70xopt=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('高斯过程回归')

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

逻辑回归(logistic regression)原理理解+matlab实现
使用梯度下降法迭代: function theta =logisticReg() % 梯度下降法寻找最合适的theta,使得代价函数J最小 options=optimset('GradObj','on ...
PID控制参数整定（调节方法）原理+图示+MATLAB调试
序首先最重要的是了解每个参数调节了系统响应的那些属性,通过观察响应从而调节参数改变属性. PID的作用概述: 1.P产生响应速度和力度,过小响应慢,过大会产生振荡,是I和D的基础. 2.I在有系统误 ...
matlab里的pid参数调节,PID控制参数整定（调节方法）原理+图示+MATLAB调试
序首先最重要的是了解每个参数调节了系统响应的那些属性,通过观察响应从而调节参数改变属性. PID的作用概述: 1.P产生响应速度和力度,过小响应慢,过大会产生振荡,是I和D的基础. 2.I在有系统误 ...
脑电图（EEG）信号去噪方法简述
前言脑电图作为目前研究最为广泛的认知大脑的方式之一,其无创性.便携性.廉价等优点都表明该方式具有巨大的发展空间.但是由于颅骨和头皮对于电信号的传输影响,从头皮采集的电信号往往混杂着非常多的噪声,并且 ...
基追踪算法 matlab,什么是基追踪算法?基于改进基追踪方法的信号去噪 - 全文
什么是基追踪算法基追踪(basis pursuit)算法是一种用来求解未知参量L1范数最小化的等式约束问题的算法. 基追踪是通常在信号处理中使用的一种对已知系数稀疏化的手段.将优化问题中的L0范数转 ...
【信号去噪】基于小波阈值实现心电信号去噪附matlab代码
1 简介由于外界环境的干扰,导致在实际信号的采集过程中无法避免地引入一些随机噪声,从而影响下一步的信号处理,所以如何对含噪信号进行去噪处理,提取出对研究有用的信号,成为信号领域的一个重要研究课题.小 ...
Python Qt GUI设计：信号与槽的使用方法（基础篇—7）
目录 1.信号与槽的概念 2.信号与槽的基础函数 2.1.创建信号函数 2.2.连接信号函数 2.3.断开信号函数 2.4.发射信号函数 3.信号和槽的使用方法 3.1.内置信号与槽的使用 3.2.自 ...
基于matlab的信号去噪,基于matlab的信号去噪研究.doc
基于matlab的信号去噪研究.doc 学生毕业设计(论文) 课题名称基于 MATLAB 的信号去噪研究姓名陈文学号 051220206 院系物电系专业电子信息工程 ...
基于MATLAB的小波收缩法信号去噪
信号在采集.传输和获取的过程中难免会受到各种噪声的干扰,这些噪声也极大地影响着人们从信号中获取有用的信息,因此对信号去噪的研究显得十分重要.小波去噪是近年来信号与信息处理领域的热门技术,利用小波算法对 ...
心电信号去噪(part5)--小波及其他技术方面问题
本系列上一篇(EMD)指路:https://blog.csdn.net/m0_37422217/article/details/90922047 小波变换我觉得这篇博客写的很好,至少我觉得,所以... ...

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