Loại bỏ tiếng ồn từ trăn tín hiệu

nội dung

• Nhập thư viện, tạo tín hiệu và thêm tiếng ồn
• Thực hiện biến đổi Fourier nhanh
• Lọc ra tiếng ồn
• Trực quan hóa kết quả
• Khử nhiễu dữ liệu thực sử dụng ngưỡng công suất
• Bộ lọc dựa trên obspy
• kết luận
• Người giới thiệu

Phân tích Fourier dựa trên ý tưởng rằng bất kỳ chuỗi thời gian nào cũng có thể được phân tách thành tổng tích phân của các sóng hài có tần số khác nhau. Do đó, về mặt lý thuyết, chúng ta có thể sử dụng một số sóng hài để tạo ra bất kỳ tín hiệu nào.

Chuỗi Fourier cho một hàm thời gian tùy ý \[f[t]\] được xác định trong khoảng [[-T/2 < t < T/2 ]] là

\begin{phương trình}

f[t] = a_0 + \sum_{n=1}^{\infty} a_n cos[\frac{2n\pi t}{T}] + \sum_{n=1}^{\infty} b_n sin[

\end{phương trình}

Trong phương trình trên, chúng ta có thể thấy rằng \[sin[\frac{2n\pi t}{T}]\] và \[cos[\frac{2n\pi t}{T}]\] tuần hoàn với . Ở đây, các giá trị lớn hơn của \[n\] tương ứng với khoảng thời gian ngắn hơn hoặc tần số cao hơn

Trong bài đăng này, chúng tôi sẽ sử dụng phân tích Fourier để lọc với giả định rằng nhiễu chồng lấp các tín hiệu trong miền thời gian nhưng không quá chồng chéo trong miền tần số

bài viết tương tự

Nhập thư viện, tạo tín hiệu và thêm tiếng ồn

import pandas as pd
import os, sys
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [10,6]
plt.rcParams.update[{'font.size': 18}]
plt.style.use['seaborn']

## Create synthetic signal
dt = 0.001
t = np.arange[0, 1, dt]
signal = np.sin[2*np.pi*50*t] + np.sin[2*np.pi*120*t] #composite signal
signal_clean = signal #copy for later comparison
signal = signal + 2.5 * np.random.randn[len[t]]
minsignal, maxsignal = signal.min[], signal.max[]

Chúng tôi đã tạo tín hiệu của mình bằng cách tính tổng hai hàm sin có tần số khác nhau [50Hz và 120Hz]. Sau đó, chúng tôi tạo ra một mảng nhiễu ngẫu nhiên và xếp chồng nhiễu đó lên tín hiệu

Thực hiện biến đổi Fourier nhanh

## Compute Fourier Transform
n = len[t]
fhat = np.fft.fft[signal, n] #computes the fft
psd = fhat * np.conj[fhat]/n
freq = [1/[dt*n]] * np.arange[n] #frequency array
idxs_half = np.arange[1, np.floor[n/2], dtype=np.int32] #first half index

Hàm fft.fft của Numpy trả về Biến đổi Fourier rời rạc một chiều bằng thuật toán Biến đổi Fourier nhanh [FFT] hiệu quả. Đầu ra của hàm rất phức tạp và chúng tôi đã nhân nó với liên hợp của nó để thu được phổ công suất của tín hiệu nhiễu. Chúng tôi đã tạo mảng tần số bằng cách sử dụng khoảng thời gian lấy mẫu [dt] và số lượng mẫu [n]

Lọc ra tiếng ồn

Trong biểu đồ trên, chúng ta có thể thấy rằng hai tần số từ tín hiệu ban đầu của chúng ta nổi bật. Bây giờ, chúng ta có thể tạo một bộ lọc có thể loại bỏ tất cả các tần số có biên độ nhỏ hơn ngưỡng của chúng ta

## Filter out noise
threshold = 100
psd_idxs = psd > threshold #array of 0 and 1
psd_clean = psd * psd_idxs #zero out all the unnecessary powers
fhat_clean = psd_idxs * fhat #used to retrieve the signal

signal_filtered = np.fft.ifft[fhat_clean] #inverse fourier transform

Trực quan hóa kết quả

## Visualization
fig, ax = plt.subplots[4,1]
ax[0].plot[t, signal, color='b', lw=0.5, label='Noisy Signal']
ax[0].plot[t, signal_clean, color='r', lw=1, label='Clean Signal']
ax[0].set_ylim[[minsignal, maxsignal]]
ax[0].set_xlabel['t axis']
ax[0].set_ylabel['Vals']
ax[0].legend[]

ax[1].plot[freq[idxs_half], np.abs[psd[idxs_half]], color='b', lw=0.5, label='PSD noisy']
ax[1].set_xlabel['Frequencies in Hz']
ax[1].set_ylabel['Amplitude']
ax[1].legend[]

ax[2].plot[freq[idxs_half], np.abs[psd_clean[idxs_half]], color='r', lw=1, label='PSD clean']
ax[2].set_xlabel['Frequencies in Hz']
ax[2].set_ylabel['Amplitude']
ax[2].legend[]

ax[3].plot[t, signal_filtered, color='r', lw=1, label='Clean Signal Retrieved']
ax[3].set_ylim[[minsignal, maxsignal]]
ax[3].set_xlabel['t axis']
ax[3].set_ylabel['Vals']
ax[3].legend[]

plt.subplots_adjust[hspace=0.4]
plt.savefig['signal-analysis.png', bbox_inches='tight', dpi=300]

Khử nhiễu dữ liệu thực sử dụng ngưỡng công suất

Tôi có bản ghi dữ liệu gia tốc kế bằng PhidgetSpatial Precision 0/0/3 Độ phân giải cao. Tôi đã chuyển đổi nó thành định dạng Miniseed để phân tích dễ dàng

# -*- coding: utf-8 -*-
# ======================================================================================================================================================
"""
Created on Thu Apr 29 12:41:26 2021

@author: Utpal Kumar [IES, Academia Sinica]
"""
# ======================================================================================================================================================

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [10,6]
plt.rcParams.update[{'font.size': 18}]
plt.style.use['seaborn']
from obspy import read
from obspy.core import UTCDateTime


otime = UTCDateTime['2021-04-18T22:14:37'] #eq origin

filenameZ = 'TW-RCEC7A-BNZ.mseed'

stZ = read[filenameZ]
streams = [stZ.copy[]]
traces = []
for st in streams:
    tr = st[0].trim[otime, otime+120]
    traces.append[tr]
    
delta = stZ[0].stats.delta
starttime = np.datetime64[stZ[0].stats.starttime]
endtime = np.datetime64[stZ[0].stats.endtime]
signalZ = traces[0].data/10**6
minsignal, maxsignal = signalZ.min[], signalZ.max[]

t = traces[0].times["utcdatetime"] 

## Compute Fourier Transform
n = len[t]
fhat = np.fft.fft[signalZ, n] #computes the fft
psd = fhat * np.conj[fhat]/n
freq = [1/[delta*n]] * np.arange[n] #frequency array
idxs_half = np.arange[1, np.floor[n/2], dtype=np.int32] #first half index
psd_real = np.abs[psd[idxs_half]] #amplitude for first half


## Filter out noise
sort_psd = np.sort[psd_real][::-1]
# print[len[sort_psd]]
threshold = sort_psd[300]
psd_idxs = psd > threshold #array of 0 and 1
psd_clean = psd * psd_idxs #zero out all the unnecessary powers
fhat_clean = psd_idxs * fhat #used to retrieve the signal

signal_filtered = np.fft.ifft[fhat_clean] #inverse fourier transform


## Visualization
fig, ax = plt.subplots[4,1]
ax[0].plot[t, signalZ, color='b', lw=0.5, label='Noisy Signal']
ax[0].set_xlabel['t axis']
ax[0].set_ylabel['Accn in Gal']
ax[0].legend[]

ax[1].plot[freq[idxs_half], np.abs[psd[idxs_half]], color='b', lw=0.5, label='PSD noisy']
ax[1].set_xlabel['Frequencies in Hz']
ax[1].set_ylabel['Amplitude']
ax[1].legend[]

ax[2].plot[freq[idxs_half], np.abs[psd_clean[idxs_half]], color='r', lw=1, label='PSD clean']
ax[2].set_xlabel['Frequencies in Hz']
ax[2].set_ylabel['Amplitude']
ax[2].legend[]

ax[3].plot[t, signal_filtered, color='r', lw=1, label='Clean Signal Retrieved']
ax[3].set_ylim[[minsignal, maxsignal]]
ax[3].set_xlabel['t axis']
ax[3].set_ylabel['Accn in Gal']
ax[3].legend[]

plt.subplots_adjust[hspace=0.6]
plt.savefig['real-signal-analysis.png', bbox_inches='tight', dpi=300]

Bộ lọc dựa trên obspy

Obspy làm cho nhiệm vụ của chúng tôi dễ dàng hơn nhiều bằng cách giới thiệu các chức năng bộ lọc. Ở đây, tôi đã sử dụng bộ lọc Butterworth-Bandpass. Để biết chi tiết về các loại bộ lọc khác nhau, bạn có thể xem tài liệu của nó

Trong ví dụ này, tôi đã sử dụng tần số góc thấp của băng thông là 0. 01 và tần số góc cao là 3 Hz dựa trên phổ tần số thu được ở trên

import pandas as pd
import os, sys
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [10,6]
plt.rcParams.update[{'font.size': 18}]
plt.style.use['seaborn']
from obspy import read
from obspy.core import UTCDateTime

otime = UTCDateTime['2021-04-18T22:14:37'] #eq origin

filenameZ = 'TW-RCEC7A-BNZ.mseed'

stZ = read[filenameZ]
streams = [stZ.copy[]]
traces = []
for st in streams:
    tr = st[0].trim[otime, otime+120]
    traces.append[tr]
    

signalZ = traces[0].data/10**6
minsignal, maxsignal = signalZ.min[], signalZ.max[]

t = np.arange[0, traces[0].stats.npts / traces[0].stats.sampling_rate, traces[0].stats.delta]

# Filtering with a lowpass on a copy of the original Trace
freqmin = 0.01
freqmax = 3
tr_filt = traces[0].copy[]
tr_filt.detrend["linear"]
tr_filt.taper[max_percentage=0.05, type='hann']
tr_filt.filter["bandpass", freqmin=freqmin,
                          freqmax=freqmax, corners=4, zerophase=True]
print[tr_filt.data/10**6]
signal_filtered = tr_filt.data/10**6

## Visualization
fig, ax = plt.subplots[2,1]
ax[0].plot[t, signalZ, color='b', lw=0.5, label='Noisy Signal']
ax[0].set_xlabel['t axis']
ax[0].set_ylabel['Accn in Gal']
ax[0].legend[]


ax[1].plot[t, signal_filtered, color='r', lw=1, label='Clean Signal Retrieved']
ax[1].set_xlabel['t axis']
ax[1].set_ylabel['Accn in Gal']
ax[1].set_title[f"Filtered in range {freqmin}-{freqmax} Hz"]
ax[1].legend[]

plt.subplots_adjust[hspace=0.4]
plt.savefig['real-signal-analysis.png', bbox_inches='tight', dpi=300]

kết luận

Trong bài viết này, chúng tôi chỉ sử dụng loại bộ lọc cơ bản để loại bỏ tiếng ồn. Với bộ lọc nâng cao, chúng tôi có thể kiểm soát nhiều hơn trong việc loại bỏ các tần số nhưng khái niệm tổng thể rất giống nhau. Trong bài đăng tiếp theo, chúng ta sẽ xem cách chúng ta có thể sử dụng wavelet để loại bỏ nhiễu

Người giới thiệu

1. Stein, S. , & Wysession, M. [2009]. Giới thiệu về địa chấn, động đất và cấu trúc trái đất