Python量化噪声卷积信号和傅里叶时频分析

🎯要点

1. 信号生成、绘制和推理图形结果：

   1. 连续时域数字信号：🎯正弦信号：🖊通用表达式、🖊不同相位式、🖊三相式。🎯指数信号：📌实数指数信号情景：如：油价指数增长函数、放射性衰减指数衰减函数，📌复数指数信号情景：傅里叶变换、正弦波形和余弦波形、线性时不变系统的特征函数，🖊实数指数、🖊PN结二极管的正向特性、🖊放射性衰减、🖊复数指数函数、🖊从复数指数提取正弦和余弦函数、🖊正弦信号与增长和衰减的实指数信号相乘。🎯非平稳信号、平稳信号、线性调频信号、方波、三角波、锯齿波、sinc 、脉冲信号和高斯信号。

   2. 离散时域信号：🎯脉冲函数列、🎯梳状函数的频域、🎯时域混叠、🎯频域混叠、🎯听觉混叠、🎯中平量化器传输特性、🎯中平量化器量化输入锯齿信号、🎯中平量化器的传输特性、🎯中升量化器量化输入锯齿信号、🎯语音信号量化、🎯均匀中平量化图像、🎯$\mu$​-法则压扩编码输入信号、🎯输入信号和量化信号间的误差、🎯量化误差概率密度函数、🎯正弦信号零阶保持插值、🎯一阶保持插值、🎯正弦信号的理想或 Sinc 插值、🎯零阶保持与 Sinc 插值的比较。🎯单位样本序列生成方式：🖊使用逻辑运算、🖊使用Python库中内置函数、🖊单元步序列、🖊单元斜坡信号、🖊实数和复数指数信号，🎯数学运算：🖊幅度缩放、🖊幅度漂移、🖊信号乘积、🖊信号叠加、🖊上/下取样、🖊时移运算、:🖊使用Python内置/无内置函数的时间反转、🖊单位样本卷积、🖊信号与移位单位样本信号的卷积、🖊卷积特征演示：交换、结合和分布、🖊低通滤波系数的方波卷积、🖊高通滤波系数的方波卷积、🎯信号相似性：🖊正弦波和余弦波的自相关和互相关、🖊正弦波与其自身和噪声信号的自相关、🖊使用自相关的延迟估计。

   3. 离散时域系统：🎯解差分方程、🎯确定状态空间转换关系、🎯计算系统的脉冲和阶跃响应、🖊计算阶跃响应、🖊从阶跃响应获取脉冲响应、🎯绘制：🖊非递归系统的零点、🖊非递归系统的幅度和相位响应

   4. 数控信号处理应用：🎯语音识别、🎯电能传输质量、🎯潘-汤普金斯算法心电图检测、🎯混合模式信号链中噪声分析。

🍇Python傅立叶分析加速度计信号去噪

傅里叶分析的思想是任何时间序列都可以分解为不同频率的谐波的积分和。因此，理论上，我们可以利用多个谐波来生成任何信号。在区间 （$-T/2<t<T/2$） 上定义的任意时间函数 f ( t ) 的傅立叶级数
$$f(t)=a_0+\sum _{n=1}^{\infty }{a}_{n}cos(\frac{2n\pi t}{T})+\sum _{n=1}^{\infty }{b}_{n}sin(\frac{2n\pi t}{T})$$
在上式中，我们可以看到 $sin(\frac{2n\pi t}{T})$ 和 $cos(\frac{2n\pi t}{T})$ 是周期性的，周期为 T。

我们将使用傅立叶分析进行滤波，假设噪声在时域中与信号重叠，但在频域中不那么重叠。

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()

我们通过对不同频率（50Hz 和 120Hz）的两个正弦函数求和来创建信号。然后我们创建了一组随机噪声并将该噪声叠加到信号上。

快速傅里叶变换

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) 创建频率数组。

在上图中，我们可以看到原始信号的两个频率很突出。现在，我们可以创建一个滤波器，可以删除幅度小于阈值的所有频率。

## 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

## 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)

功率阈值对真实数据进行去噪

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)

参阅一：计算思维

参阅二：亚图跨际