如何用Python区分信号的下降沿？
我正在处理一些信号，以计算一个弹簧-质量系统（地震仪）的自由振荡周期和阻尼比。我使用Python作为主要的处理程序。我需要做的就是导入这个信号，解析信号并找到下降沿，然后返回一个包含每次振荡顶峰的列表，这样我就可以计算阻尼比。一旦我确定了数据集中振荡的位置，计算自由周期就相对简单了。
我现在主要卡在如何解析这个列表，识别下降沿，然后捕捉每个Z0..Zn元素上。振荡频率可以通过快速傅里叶变换（FFT）相对容易地计算出来，一旦我知道了下降沿的位置，但如果我处理整个文件，系统在释放前的能量可能会导致算法输出一个超低频率，这个频率代表的是近直流偏移，而不是实际的振荡。（在高阻尼比的情况下，这尤其成问题，因为可能只有四五次可测的反弹）。
有没有人有什么想法可以帮我解决这个问题？现在，下面的代码使用了截图中信号的任意赋值。但是我需要代码来计算这些值。此外，我还没有确定如何创建我的峰值列表以计算阻尼比h。对于解决这个问题，如果你能提供一些想法，我将非常感激。由于我在Stackoverflow上的声誉很小，我在以下链接中包含了我的信号的示例截图： （希望这个链接能用！） Tycho的示例信号 -->
##########################################################
import os, sys, csv
from scipy import signal
from scipy.integrate import simps
import pylab as plt
import numpy as np
import scipy as sp
#
# code goes here that imports the data from the csv file into lists.
# One of those lists is called laser, a time-history list in counts from an analog-digital-converter
# sample rate is 130.28 samples / second.
#
#
# Find the period of the observed signal
#
delta = 0.00767 # Calculated elsewhere in code, represents 130.28 samples/sec
# laser is a list with about 20,000 elements representing time-history data from a laser position sensor
# The release of the mass and system response occurs starting at sample number 2400 in this particular instance.
sense = signal.detrend(laser[2400:(2400+8192)]) # Remove the mean of the signal
N = len(sense)
W = np.fft.fft(sense)
freq = np.fft.fftfreq(len(sense),delta) # First value represents the number of samples and delta is the sample rate
#
# Take the sample with the largest amplitude as our center frequency.
# This only works if the signal is heavily sinusoidal and stationary
# in nature, like our calibration data.
#
idx = np.where(abs(W)==max(np.abs(W)))[0][-1]
Frequency = abs(freq[idx]) # Frequency in Hz
period = 1/(Frequency*delta) # represents the number of samples for one cycle of the test signal.
#
# create an axis representing time.
#
dt = [] # Create an x axis that represents elapsed time in seconds. delta = seconds per sample, i represents sample count
for i in range(0,len(sensor)):
dt.append(i*delta)
#
# At this point, we know the frequency interval, the delta, and we have the arrays
# for signal and laser. We can now discriminate out the peaks of each rebound and use them to process the damping ratio of either the 'undamped' system or the damping ratio of the 'electrically damped' system.
#
print 'Frequency calcuated to ',Frequency,' Hz.'
1 个回答
这里有一个有点特别的想法，我觉得这个方法可能很有效，而且不需要太多的猜测和经验法则。你手里的数据质量很高，而且符合一个已知的曲线，这对我们来说是个好消息。在这里，我假设你曲线的“好部分”是这样的：
V = a * exp(-γ * t) * cos(2 * π * f * t + φ) + V0 # [Eq1]
V: 电压t: 时间γ: 阻尼常数f: 频率a: 起始幅度φ: 起始相位V0: 直流偏移
算法概述
去掉偏移
首先，进行数值求导。因为数据质量相当高，所以噪声对结果的影响应该不大。
电压的导数 V_deriv 和原始数据的形式是一样的：频率和阻尼常数相同，但相位 ψ 和幅度 b 不同，
V_deriv = b * exp(-γ * t) * cos(2 * π * f * t + ψ) # [Eq2]
好处是，这样可以自动去掉你的直流偏移。
替代方案：这一步并不是绝对必要的——偏移对曲线拟合来说只是一个相对小的复杂因素，因为你总是可以通过对整个曲线取平均来提供一个好的偏移估计。这样做的权衡是，如果不使用导数，你会得到更好的信噪比。
初步猜测
现在考虑导数的曲线（如果你跳过了上一步，就看原始曲线）。从最右边的数据点开始，然后向左跟随曲线，尽可能多地找到振荡，直到你找到一个振荡，其幅度大于某个幅度阈值 A。你想要找到一个包含一个振荡的曲线部分，并且信噪比良好。
如何确定幅度阈值比较难说。这取决于你的传感器有多好。我建议把这个作为一个参数，方便后续调整。
曲线拟合
现在你已经捕捉到一个振荡，可以做很多简单的事情：估计频率 f 和阻尼常数 γ。利用这些初步估计，你可以对这个振荡右侧的所有数据进行非线性曲线拟合（包括这个振荡本身）。
如果你使用导数，拟合的函数是[Eq2]，如果使用原始曲线则是[Eq1]（不过这两个函数其实是一样的，所以这个区别不大）。
为什么需要这些初步估计呢？对于衰减波的非线性曲线拟合，给出好的参数初始猜测是很关键的，尤其是频率和阻尼常数。其他参数相对不那么重要（至少根据我的经验），但如果你需要获取其他参数，可以这样做：
- 如果你从最大值开始，相位是
0；如果从最小值开始，相位是π。 - 你也可以猜测起始幅度，尽管拟合算法通常即使你设置为
1也能很好地工作。
找到边界
此时曲线拟合应该非常好——你可以通过检查你的曲线和拟合之间的差异来判断，因为这个差异非常小。现在的关键是尝试向左扩展曲线的范围。
如果你保持在正弦波区域，曲线拟合应该会保持相当好（如何判断“好”需要一些实验）。然而，一旦你进入曲线的平坦区域，误差会开始显著增加，参数也会开始偏离。这可以用来确定“好”数据的结束位置。
不过你不必逐点进行，这样效率很低。可以使用二分查找，这在这里应该效果很好（可能是它的“单边”变体）。
总结
你不必完全按照这个流程，但基本的意思是，你可以对从最右边开始的一小部分数据进行一些分析，然后逐渐扩展时间范围，直到你发现误差开始增大而不是减小。
此外，你还可以结合不同的经验法则，看看它们是否一致。如果不一致，那可能是一个比较可疑的数据点，需要一些手动干预。
需要注意的是，我上面提到的算法的一个特别优点是，你会在这个过程中得到你想要的结果（阻尼常数和频率），同时还会有不确定性估计。
我没有提到大部分复杂的数学和算法细节，以便提供一个总体概述，但如果需要，我可以提供更多细节。