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

import numpy as npdef 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

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 npdef 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 signaldef moving_average_sg(data, window_size, order):        # SG滤波    sg_data = signal.savgol_filter(data, window_size, order)        return sg_data

import numpy as npfrom numpy.linalg import svddef 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

import matplotlib.pyplot as pltimport 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)