在处理时域信号时,常常会遇到数据含有噪音的情况,这在坐标图上的表现就是图像是有很多微x小的的“突刺”;对于一些精度要求偏高的情形,会影响后续信号特征的提取以及频域信号的转换。
我们的任务是尽量让图像尽可能地光滑,而它的一些阶跃点我们仍然保留。基于这样的思想,可以用Tikhonov正则化方法来实现这一任务。
对于一般的拟合任务,我们假定x为去噪信号,y为原始信号,那么我们既要希望去噪信号与原始信号之间的差异尽可能地小,同时去噪信号x还要尽可能地光滑一些。因此我们让正则化项为x的二阶导数,这里用二阶的梯度算符表示,而拟合项为x与y的差(K取单位矩阵)。这样一来,通过选取合适的正则化参数α,我们就可以兼顾去噪信号x的光滑性,和基本特征的保留程度。
如此一来,就有:
事实上,当我们选定了合适的正则化参数以后,泛函就是已经确定了的,这样的话,我们只需要找到使得该函数最小化的x,这就是得到的去噪信号x,具体的求解
的最小值一般有两种方法,一种是令
,将
用二阶中心差分近似代替(步长的系数由正则化系数α决定),这样得到的三对角矩阵记为L,即
.
又因为,
,
于是,用最小二乘法解出这个方程中的x,就得到了去噪的信号x。实际求解过程中,使用默认的高斯消元法求解会消耗MATLAB较多的内存和算力,并且求解速度较慢于第二种做法;所以如果要使用本方法,需要配合雅可比迭代法,或者高斯-赛德尔方法快速实现矩阵运算。
而另一种方式是用梯度下降(Gradient Descent)算法, 具体的原理这里略去。总而言之,使用这种算法,是不断地迭代去噪信号x,直到迭代后的x最终满足一定的条件,此时即认为得到的x为去噪信号。本文使用的matlab代码即是使用方法二,并在此基础上,使用了变种Nestrov动量法来加速梯度下降的过程。即在梯度下降的过程中,受前一步的影响越来越大,改良算法可以更快地达到收敛。
具体的Nestrov算法是这样的:
,
以下是Matlab代码
function [dinoised_signal,raw_signal] = Tikhonov_Nestrov(signal,eps,h,alpha)
%%INPUT
%signal 去噪信号向量
%eps 精度,当精度达到某要求时,返回所得到的去噪信号
%h 改良的梯度优化算法的步长值,可以设置,以更好地更快地实现去噪
%alpha 正则化系数,显示比重
%%OUTPUT
if size(signal,1) == 1
signal = signal';%转置为列向量
end
n = size(signal,1);
L = zeros(n-2,n);
K = eye(n);
for idx = 1:n-2
L(idx,idx) = 1;
L(idx,idx+1) = -2;
L(idx,idx+2) = 1;
end
f = @(x).5*norm(K*x-signal)+.5*alpha*norm(L*x);
df = @(x)K'*(K*x-signal)+alpha*(L'*L)*x;
iter = 1;
x0 = signal;
grad = df(x0);
x = x0;
U = zeros(n,2);
while (any(grad>eps))
U(:,1) = x - h*grad;
x = U(:,1) + (iter-1)/(iter+2)*(U(:,1)-U(:,2));
U(:,2) = U(:,1);
% 上面三行是改进的Nestrov动量法
% x = x-h*grad;%如果把上面三行代码换成这一行,就是最基本的梯度下降法
grad = df(x);
iter = iter+1;
% if iter>=1000
% x = x0;
% break;
% end
end
%这里的最终的x即为去噪信号,而x0是原始信号,可以用绘图的方式比较两者的差异
dinoised_signal = x;raw_signal = x0;
end
具体的效果可以参考下面一张图,这是笔者在一次数模论文中的插图,红色是原始信号,黑色是去噪信号。选取不同的正则化系数有不同的效果,一般为0.001,0.01,0.1,不会太大
具体的滤波降噪算法有很多种,还有比较常见的就是中值滤波算法,卡尔曼滤波算法。对于不同的实际问题,选用不同的滤波算法收获的效果是完全不同的,其他的滤波算法,未来有空我也会继续分享。以上代码都是可以轻松移植并使用的,如果本文对你有帮助的话,请点一个免费的赞,您的支持是我继续更新的动力!