信号处理基础之噪声与降噪(二) | 平滑降噪、SVD降噪及其python代码实现

信号侦探社

于 2024-05-11 15:49:03 发布

阅读量1.7k

点赞数 43

分类专栏：信号处理文章标签：信号处理 python 开发语言

本文链接：https://blog.csdn.net/weixin_43394610/article/details/138717368

版权

信号处理专栏收录该内容

4 篇文章 2 订阅

订阅专栏


import numpy as np
def denoise_evaluate(signal_pure, noise, signal_denoised):
    value = {}    # SNR计算
    p_signal = np.sum(signal_pure**2) / len(signal_pure)    
    p_noise = np.sum(noise**2) / len(noise)    
    SNR = 10 * np.log10(p_signal / p_noise)    
    value['SNR'] = SNR    
    # NCC计算    
    NCC = np.correlate(signal_pure, signal_denoised) / (np.linalg.norm(signal_pure) * np.linalg.norm(signal_denoised))    value['NCC'] = NCC[0]    
    # MSE计算    
    MSE = np.mean((signal_pure - signal_denoised) ** 2)    
    value['MSE'] = MSE    
    # RMSE计算    
    RMSE = np.sqrt(MSE)    
    value['RMSE'] = RMSE    
    return value

注：上述四个指标的计算均依赖于”纯净信号“，对于只知道真是含噪信号，不知道纯净信号的情况，无法进行计算，即对于工程实际中的随机信号，无法直接用上述四个指标进行衡量。

2 平滑降噪

2.1 滑动平均的基本原理

滑动平均（Moving Average）是一种时域滤波方法，其基本原理是设定一个滑动窗口，该滑动窗口沿着原始数据时序方向移动，每次移动时计算当前窗口的平均值作为滤波值，最终得到滑动平均序列，其过程可表示为下图。

记滑动窗口长度为N（N为奇数），则滑动平均可表示为：

$y(n)=\cfrac{1}{N}(x(n-\cfrac{N-1}{2})+\cdots+x(n)+\cdots+x(n+\cfrac{N-1}{2}))=\sum_{k=-\cfrac{N-1}{2}}^{\cfrac{N-1}{2}}x(n+k)$

滑动窗口为对称窗口，防止出现相位偏差，窗口长度一般为奇数。特别地，滑动平均的滤波效果取决于滑动窗的长度，一般长度越大平滑效果越好，但窗口选择过大可能出现过平滑，湮没有效信号，且滑动平均滤波对边缘数据的处理效果不佳。

注：也可通过卷积的方式实现滑动滤波，此处不展开讨论

2.2 SG滤波基本原理

SG(Savitzky-Golay)滤波是一种广泛使用的数据平滑滤波降噪方法，其基于曲线局部特征进行多项式拟合，应用最小二乘法确定加权系数进行移动窗口加权平均，重构的数据能够较好保留局部特征，不受时间及空间尺度的影响。SG滤波的基本公式可表示为：

$y_j^*=\cfrac{\sum_{i=-m}^{i=m}C_iY_{j+i}}{N}$

具体地，记滤波窗口长度为 N（ N为奇数），测量点记为 $\bold{\mathit{x}}=(x_1,x_2,\cdots,x_n)$ ，采用k-1次多项式对窗口内的数据点进行拟合，记拟合常数为 $\bold{\mathit{A}}=(a_0,a_1,\cdots,a_{k-1})$ ，则有滤波结果

$y=\bold{\mathit{A}}\cdot[1,\bold{\mathit{x}},\bold{\mathit{x}}^2,\cdots,\bold{\mathit{x}}^{k-1}]^T$

不难发现，在使用SG滤波进行降噪时，窗口长度N，多项式拟合阶次k对降噪结果起着至关重要的作用。窗口长度越长，单次参与拟合的数据量越多，多项式拟合阶次越高，单次拟合的结果越接近原始信号，保留的细节越多，但是过高的阶次可能出现过拟合，产生新的噪声。

注：限于篇幅，不对完整过程进行推导，详细过程请参考文献[3]。

2.3 python实现

本文给出了3种方法，常规方法滤波、卷积平滑滤波、SG滤波。

def moving_average(data, window_size):
    # 常规方法进行平滑滤波    
    smoothed_data = []    
    for i in range(len(data)):        
    start = max(0, i - window_size // 2)        
    end = min(len(data), i + window_size // 2 + 1)        
    window = data[start:end]        
    smoothed_data.append(sum(window) / len(window))    
    return smoothed_data

import numpy as np
def moving_average_conv(data, window_size):    
# 卷积平滑滤波    
    window = np.ones(window_size) / window_size    
    smoothed_data = np.convolve(data, window, mode='same')    
    return smoothed_data


import numpy as npfrom scipy 
import signal
def moving_average_sg(data, window_size, order):    
    # SG滤波
    sg_data = signal.savgol_filter(data, window_size, order)    
    return sg_data

3 SVD降噪

3.1 SVD降噪的基本原理

奇异值分解(Singular Value Decomposition, SVD)能够对矩阵进行分解，得到代表矩阵最本质变化的矩阵元素。作为信号处理的重点研究方向之一，其主要功能是去除信号中的随机噪声，降噪后的信号不存在时间延迟且相移较小。

记一个含噪声的时间序列 $y(t)$ ， $y(t)$ 表示为理想信号 $x(t)$ 和噪声信号 $r(t)$ 的和，表示为

$y(t)=x(t)+r(t)$

将 $y(t)=[y_1,y_2,\cdots,y_N]$ 构造成Hankel矩阵（含噪矩阵A），表示为：

$\bold{\mathit{A}}=\begin{bmatrix}y_1&y_2&\cdots&y_n\\y_2&y_3&\cdots&y_{n+1}\\\vdots&\vdots&\vdots&\vdots\\y_m&y_{m+1}&\cdots&y_{m+n-1}\end{bmatrix}$

式中，当N为偶数时， $n=N/2$ ，N为奇数时， $n=(N+1)/2$ ， $N=m+n-1$ 。

对含噪矩阵A进行奇异值分解可得：

$\bold{\mathit{A}}=\bold{\mathit{USV^T}}$

$\bold{\mathit{S}}=diag(\sigma_1,\sigma_2,\cdots,\sigma_q,\bold{\mathit{O}})\,\,q=min(m,n)$

其中 $\sigma_i$ 为矩阵A的奇异值， $\sigma_1>\sigma_2>\cdots>\sigma_q$ ，O为零矩阵，U、V为正交矩阵。

选取有效奇异值个数p，保留主对角矩阵S的前p个较大的奇异值，将其余奇异值全部置为0，得到新的对角矩阵：

$\bold{\mathit{S'}}=diag(\sigma_1,\sigma_2,\cdots,\sigma_p,\bold{\mathit{O}})$

于是可以通过左右正交矩阵U和V得到新的信号矩阵 $A'$ ，将重构的矩阵 $A'$ 恢复成一维信号，此时便完成了SVD降噪处理。

3.2 python实现

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft
from numpy.linalg import svd

def denoised_with_svd(data, nlevel=8):
    """
    :param data: 需要降噪的原始数据，1D-array
    :param nlevel: 阶次
    :return:重构后的信号
    """
    N = len(data)
    A = np.lib.stride_tricks.sliding_window_view(data, (N//2,))
    U, S, Vh = svd(A)

    # 重构信号
    X = np.zeros_like(A)
    for i in range(nlevel):
        X += Vh[i, :] * S[i] * U[:, i:i+1]
    X = X.T
    data_recon = np.zeros(N)
    for i in range(N):
        a = 0
        m = 0
        for j1 in range(N//2):
            for j2 in range(N//2+1):
                if i == j1+j2:
                    a = a + X[j1, j2]  # 把矩阵重构回一维信号
                    m = m + 1
    
        if m != 0:
            data_recon[i] = a / m
    return data_recon

4 算例测试

4.1 测试用例

import matplotlib.pyplot as plt
import matplotlib
# 测试用例
n = 500  # 生成500个点的信号
t = np.linspace(0, 30*np.pi, n, endpoint=False)
s = np.cos(0.1*np.pi*t) + np.sin(0.3*np.pi*t) + np.cos(0.5*np.pi*t) + np.sin(0.7*np.pi*t)  # 原始信号
r = np.random.randn(n)
y = s + r  # 加噪声
data_smooth_avg = moving_average(y, 5)
data_smooth_conv = moving_average_conv(y, 5)
data_smooth_sg = moving_average_sg(y, 5, 3)
data_svd = denoised_with_svd(y)
fig = plt.figure(figsize=(16, 9))
plt.subplots_adjust(hspace=0.5)
font = {'family': 'Times New Roman', 'size': 16, 'weight': 'normal',}
matplotlib.rc('font', **font)
plt.subplot(6, 1, 1)
plt.plot(s, linewidth=1.5, color='b')
plt.title('原始模拟信号', fontname='microsoft yahei', fontsize=16)
plt.subplot(6, 1, 2)
plt.plot(y, linewidth=1.5, color='r')
plt.title('原始模拟信号+高斯噪声', fontname='microsoft yahei', fontsize=16)
plt.subplot(6, 1, 3)
plt.plot(data_smooth, linewidth=1.5, color='g')
plt.title('降噪信号-均值法', fontname='microsoft yahei', fontsize=16)
plt.subplot(6, 1, 4)
plt.plot(data_smooth, linewidth=1.5, color='g')
plt.title('降噪信号-卷积法', fontname='microsoft yahei', fontsize=16)
plt.subplot(6, 1, 5)
plt.plot(data_smooth, linewidth=1.5, color='g')
plt.title('降噪信号-SG滤波', fontname='microsoft yahei', fontsize=16)
plt.subplot(6, 1, 6)
plt.plot(data_svd, linewidth=1.5, color='g')
plt.title('降噪信号-SVD', fontname='microsoft yahei', fontsize=16)