简单的基于小波分析的时间序列降噪方法(Python)

import numpy as np
import pywt
import matplotlib.pyplot as plt


def denoise_signal(signal, wavelet='db4', level=6, noise_std=0.1):
    try:
        # Decompose signal into wavelet coefficients
        coeffs = pywt.wavedec(signal, wavelet, level=level)


        # Threshold the coefficients using universal threshold
        threshold = noise_std * np.sqrt(2 * np.log(len(signal)))
        denoised_coeffs = [pywt.threshold(c, threshold, mode='soft') for c in coeffs]


        # Reconstruct the signal from the denoised coefficients
        denoised_signal = pywt.waverec(denoised_coeffs, wavelet)


        # Ensure the reconstructed signal has the same length as the original signal
        denoised_signal = denoised_signal[:len(signal)]


        # Calculate SNR
        noise = signal - denoised_signal
        snr = 20 * np.log10(np.linalg.norm(signal) / np.linalg.norm(noise))


        return denoised_signal, snr
    except Exception as e:
        print(f"Error during denoising: {e}")
        return None, None


def plot_signals(original_signal, noisy_signal, denoised_signal):
    plt.figure(figsize=(12, 6))
    plt.plot(original_signal, label='Original signal')
    plt.plot(noisy_signal, label='Noisy signal')
    plt.plot(denoised_signal, label='Denoised signal')
    plt.legend()
    plt.title('Signal Denoising using Wavelet Transform')
    plt.xlabel('Sample Index')
    plt.ylabel('Amplitude')
    plt.grid(True)
    plt.show()


# Generate a noisy signal
np.random.seed(0)
original_signal = np.random.randn(1000)
noise = 0.1 * np.random.randn(1000)
noisy_signal = original_signal + noise


# Denoise the signal
denoised_signal, snr = denoise_signal(noisy_signal)




plot_signals(original_signal, noisy_signal, denoised_signal)




if snr is not None:
    print('SNR:', snr)
else:
    print('Denoising failed.')

import numpy as np
import pywt


def add_noise(signal, noise_std):
    noise = np.random.normal(0, noise_std, size=len(signal))
    noisy_signal = signal + noise
    return noisy_signal


def denoise_signal(signal):
    # Decompose signal into wavelet coefficients
    coeffs = pywt.wavedec(signal, 'db4', level=6)


    # Threshold the coefficients using universal threshold
    threshold = np.sqrt(2*np.log(len(signal)))
    denoised_coeffs = [pywt.threshold(c, threshold, mode='soft') for c in coeffs]


    # Reconstruct the signal from the denoised coefficients
    denoised_signal = pywt.waverec(denoised_coeffs, 'db4')


    # Calculate the signal-to-noise ratio (SNR)
    noise = signal - denoised_signal
    snr = 20*np.log10(np.linalg.norm(signal)/np.linalg.norm(noise))


    return denoised_signal, snr


# Generate a signal
signal = np.sin(2*np.pi*5*np.linspace(0, 1, num=1000))


# Add Gaussian noise to the signal
noise_std = 0.1
noisy_signal = add_noise(signal, noise_std)


# Denoise the signal and print the SNR
denoised_signal, snr = denoise_signal(noisy_signal)
print('SNR:', snr)

SNR: 3.5626499602068638

知乎学术咨询：https://www.zhihu.com/consult/people/792359672131756032?isMe=1

担任《Mechanical System and Signal Processing》等审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。