关键词Python | 时域 | 信号 | 连续 | 离散 | 相位 | 指数 | 函数 | 复数 | 正弦 | 余弦 | 波形 | 线性 | 时不变系统 | 特征函数 | 衰减 | 调频 | 非平稳 | 平稳 | 脉冲 | 梳状函数 | 混叠 | 听觉 | 量化器 | 图像 | 样本 | 幅度 | 漂移 | 滤波 | 语音 | 电能 | 质量 | 变换 | 系数 | 噪音
傅里叶分析的思想是任何时间序列都可以分解为不同频率的谐波的积分和。因此,理论上,我们可以利用多个谐波来生成任何信号。在区间 (- T / 2< t < T / 2) 上定义的任意时间函数 f ( t ) 的傅立叶级数
Copy 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')
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()
Copy n = len(t)
fhat = np.fft.fft(signal, n)
psd = fhat * np.conj(fhat)/n
freq = (1/(dt*n)) * np.arange(n)
idxs_half = np.arange(1, np.floor(n/2), dtype=np.int32)
Numpy 的 fft.fft 函数使用高效的快速傅里叶变换 (FFT) 算法返回一维离散傅里叶变换。 该函数的输出是复数,我们将其与其共轭相乘以获得噪声信号的功率谱。 我们使用采样间隔 (dt) 和样本数 (n) 创建频率数组。
Copy ## 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
Copy ## 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)
Copy 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')
otime = UTCDateTime('2021-04-18T22:14:37')
filenameZ = 'TRCEC7A.msd'
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")
n = len(t)
fhat = np.fft.fft(signalZ, n)
psd = fhat * np.conj(fhat)/n
freq = (1/(delta*n)) * np.arange(n)
idxs_half = np.arange(1, np.floor(n/2), dtype=np.int32)
psd_real = np.abs(psd[idxs_half])
sort_psd = np.sort(psd_real)[::-1]
threshold = sort_psd[300]
psd_idxs = psd > threshold
psd_clean = psd * psd_idxs
fhat_clean = psd_idxs * fhat
signal_filtered = np.fft.ifft(fhat_clean)
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)