在傅立叶或时域中积分
Integration in Fourier or time domain
我很难理解信号的数值积分问题。基本上我有一个信号，我想集成或执行和反微分作为时间的函数（集成拾取线圈以获得磁场）。我尝试了两种不同的方法，原则上应该是一致的，但事实并非如此。我正在使用的代码如下。请注意，代码中的信号 y 之前已使用巴特沃斯滤波进行了高通滤波（类似于此处所做的http://wiki.scipy.org/Cookbook/ButterworthBandpass). The signal and time basis can be downloaded here (https://www.dropbox.com/s/fi5z38sae6j5410/trial.npz?dl=0）
import scipy as sp
from scipy import integrate
from scipy import fftpack
data = np.load('trial.npz')
y = data['arr_1'] # this is the signal
t = data['arr_0']
# integration using pfft
bI = sp.fftpack.diff(y-y.mean(),order=-1)
bI2= sp.integrate.cumtrapz(y-y.mean(),x=t)
现在这两个信号（除了可以取出的最终不同的线性趋势）是不同的，或者更好的是动态它们非常相似，振荡时间相同，但两个信号之间的系数大约为 30 ，在 bI2 比 bI 低 30 倍（大约）的意义上。顺便说一句，我已经减去两个信号的均值以确保它们是零均值信号并在 IDL 中执行积分（均具有等效的 cumsumtrapz 和在傅立叶域中）给出与 bI2 兼容的值。真的欢迎任何线索
很难知道 scipy.fftpack.diff() 在幕后做了什么。
为了尝试解决你的问题，我挖出了我前一段时间写的一个旧的频域积分函数。值得指出的是，在实践中，人们通常希望对一些参数的控制比 scipy.fftpack.diff() 给你的要多一些。例如，我的 intf() 函数的 f_lo 和 f_hi 参数允许您 band-limit 输入以排除可能有噪声的非常低或非常高的频率。在积分期间，尤其是嘈杂的低频可以 'blow-up' 并淹没信号。您可能还想在时间序列的开始和结束处使用 window 来阻止频谱泄漏。
我已经计算了 bI2 和一个结果，bI3，使用以下代码与 intf() 集成一次（为简单起见，我假设平均采样率）：
import intf
from scipy import integrate
data = np.load(path)
y = data['arr_1']
t = data['arr_0']
bI2= sp.integrate.cumtrapz(y-y.mean(),x=t)
bI3 = intf.intf(y-y.mean(), fs=500458, f_lo=1, winlen=1e-2, times=1)
我绘制了 bI2 和 bI3：
这两个时间序列具有相同的数量级，并且形状大致相同，尽管 bI2 中存在明显的分段线性趋势。我知道这并不能解释 scipy 函数中发生的事情，但至少这表明它不是频域方法的问题。
下面完整粘贴了 intf 的代码。
def intf(a, fs, f_lo=0.0, f_hi=1.0e12, times=1, winlen=1, unwin=False):
"""
Numerically integrate a time series in the frequency domain.
This function integrates a time series in the frequency domain using
'Omega Arithmetic', over a defined frequency band.
Parameters
----------
a : array_like
Input time series.
fs : int
Sampling rate (Hz) of the input time series.
f_lo : float, optional
Lower frequency bound over which integration takes place.
Defaults to 0 Hz.
f_hi : float, optional
Upper frequency bound over which integration takes place.
Defaults to the Nyquist frequency ( = fs / 2).
times : int, optional
Number of times to integrate input time series a. Can be either
0, 1 or 2. If 0 is used, function effectively applies a 'brick wall'
frequency domain filter to a.
Defaults to 1.
winlen : int, optional
Number of seconds at the beginning and end of a file to apply half a
Hanning window to. Limited to half the record length.
Defaults to 1 second.
unwin : Boolean, optional
Whether or not to remove the window applied to the input time series
from the output time series.
Returns
-------
out : complex ndarray
The zero-, single- or double-integrated acceleration time series.
Versions
----------
1.1 First development version.
Uses rfft to avoid complex return values.
Checks for even length time series; if not, end-pad with single zero.
1.2 Zero-means time series to avoid spurious errors when applying Hanning
window.
"""
a = a - a.mean() # Convert time series to zero-mean
if np.mod(a.size,2) != 0: # Check for even length time series
odd = True
a = np.append(a, 0) # If not, append zero to array
else:
odd = False
f_hi = min(fs/2, f_hi) # Upper frequency limited to Nyquist
winlen = min(a.size/2, winlen) # Limit window to half record length
ni = a.size # No. of points in data (int)
nf = float(ni) # No. of points in data (float)
fs = float(fs) # Sampling rate (Hz)
df = fs/nf # Frequency increment in FFT
stf_i = int(f_lo/df) # Index of lower frequency bound
enf_i = int(f_hi/df) # Index of upper frequency bound
window = np.ones(ni) # Create window function
es = int(winlen*fs) # No. of samples to window from ends
edge_win = np.hanning(es) # Hanning window edge
window[:es/2] = edge_win[:es/2]
window[-es/2:] = edge_win[-es/2:]
a_w = a*window
FFTspec_a = np.fft.rfft(a_w) # Calculate complex FFT of input
FFTfreq = np.fft.fftfreq(ni, d=1/fs)[:ni/2+1]
w = (2*np.pi*FFTfreq) # Omega
iw = (0+1j)*w # i*Omega
mask = np.zeros(ni/2+1) # Half-length mask for +ve freqs
mask[stf_i:enf_i] = 1.0 # Mask = 1 for desired +ve freqs
if times == 2: # Double integration
FFTspec = -FFTspec_a*w / (w+EPS)**3
elif times == 1: # Single integration
FFTspec = FFTspec_a*iw / (iw+EPS)**2
elif times == 0: # No integration
FFTspec = FFTspec_a
else:
print 'Error'
FFTspec *= mask # Select frequencies to use
out_w = np.fft.irfft(FFTspec) # Return to time domain
if unwin == True:
out = out_w*window/(window+EPS)**2 # Remove window from time series
else:
out = out_w
if odd == True: # Check for even length time series
return out[:-1] # If not, remove last entry
else:
return out
我很难理解信号的数值积分问题。基本上我有一个信号，我想集成或执行和反微分作为时间的函数（集成拾取线圈以获得磁场）。我尝试了两种不同的方法，原则上应该是一致的，但事实并非如此。我正在使用的代码如下。请注意，代码中的信号 y 之前已使用巴特沃斯滤波进行了高通滤波（类似于此处所做的http://wiki.scipy.org/Cookbook/ButterworthBandpass). The signal and time basis can be downloaded here (https://www.dropbox.com/s/fi5z38sae6j5410/trial.npz?dl=0）
import scipy as sp
from scipy import integrate
from scipy import fftpack
data = np.load('trial.npz')
y = data['arr_1'] # this is the signal
t = data['arr_0']
# integration using pfft
bI = sp.fftpack.diff(y-y.mean(),order=-1)
bI2= sp.integrate.cumtrapz(y-y.mean(),x=t)
现在这两个信号（除了可以取出的最终不同的线性趋势）是不同的，或者更好的是动态它们非常相似，振荡时间相同，但两个信号之间的系数大约为 30 ，在 bI2 比 bI 低 30 倍（大约）的意义上。顺便说一句，我已经减去两个信号的均值以确保它们是零均值信号并在 IDL 中执行积分（均具有等效的 cumsumtrapz 和在傅立叶域中）给出与 bI2 兼容的值。真的欢迎任何线索
很难知道 scipy.fftpack.diff() 在幕后做了什么。
为了尝试解决你的问题，我挖出了我前一段时间写的一个旧的频域积分函数。值得指出的是，在实践中，人们通常希望对一些参数的控制比 scipy.fftpack.diff() 给你的要多一些。例如，我的 intf() 函数的 f_lo 和 f_hi 参数允许您 band-limit 输入以排除可能有噪声的非常低或非常高的频率。在积分期间，尤其是嘈杂的低频可以 'blow-up' 并淹没信号。您可能还想在时间序列的开始和结束处使用 window 来阻止频谱泄漏。
我已经计算了 bI2 和一个结果，bI3，使用以下代码与 intf() 集成一次（为简单起见，我假设平均采样率）：
import intf
from scipy import integrate
data = np.load(path)
y = data['arr_1']
t = data['arr_0']
bI2= sp.integrate.cumtrapz(y-y.mean(),x=t)
bI3 = intf.intf(y-y.mean(), fs=500458, f_lo=1, winlen=1e-2, times=1)
我绘制了 bI2 和 bI3：
这两个时间序列具有相同的数量级，并且形状大致相同，尽管 bI2 中存在明显的分段线性趋势。我知道这并不能解释 scipy 函数中发生的事情，但至少这表明它不是频域方法的问题。
下面完整粘贴了 intf 的代码。
def intf(a, fs, f_lo=0.0, f_hi=1.0e12, times=1, winlen=1, unwin=False):
"""
Numerically integrate a time series in the frequency domain.
This function integrates a time series in the frequency domain using
'Omega Arithmetic', over a defined frequency band.
Parameters
----------
a : array_like
Input time series.
fs : int
Sampling rate (Hz) of the input time series.
f_lo : float, optional
Lower frequency bound over which integration takes place.
Defaults to 0 Hz.
f_hi : float, optional
Upper frequency bound over which integration takes place.
Defaults to the Nyquist frequency ( = fs / 2).
times : int, optional
Number of times to integrate input time series a. Can be either
0, 1 or 2. If 0 is used, function effectively applies a 'brick wall'
frequency domain filter to a.
Defaults to 1.
winlen : int, optional
Number of seconds at the beginning and end of a file to apply half a
Hanning window to. Limited to half the record length.
Defaults to 1 second.
unwin : Boolean, optional
Whether or not to remove the window applied to the input time series
from the output time series.
Returns
-------
out : complex ndarray
The zero-, single- or double-integrated acceleration time series.
Versions
----------
1.1 First development version.
Uses rfft to avoid complex return values.
Checks for even length time series; if not, end-pad with single zero.
1.2 Zero-means time series to avoid spurious errors when applying Hanning
window.
"""
a = a - a.mean() # Convert time series to zero-mean
if np.mod(a.size,2) != 0: # Check for even length time series
odd = True
a = np.append(a, 0) # If not, append zero to array
else:
odd = False
f_hi = min(fs/2, f_hi) # Upper frequency limited to Nyquist
winlen = min(a.size/2, winlen) # Limit window to half record length
ni = a.size # No. of points in data (int)
nf = float(ni) # No. of points in data (float)
fs = float(fs) # Sampling rate (Hz)
df = fs/nf # Frequency increment in FFT
stf_i = int(f_lo/df) # Index of lower frequency bound
enf_i = int(f_hi/df) # Index of upper frequency bound
window = np.ones(ni) # Create window function
es = int(winlen*fs) # No. of samples to window from ends
edge_win = np.hanning(es) # Hanning window edge
window[:es/2] = edge_win[:es/2]
window[-es/2:] = edge_win[-es/2:]
a_w = a*window
FFTspec_a = np.fft.rfft(a_w) # Calculate complex FFT of input
FFTfreq = np.fft.fftfreq(ni, d=1/fs)[:ni/2+1]
w = (2*np.pi*FFTfreq) # Omega
iw = (0+1j)*w # i*Omega
mask = np.zeros(ni/2+1) # Half-length mask for +ve freqs
mask[stf_i:enf_i] = 1.0 # Mask = 1 for desired +ve freqs
if times == 2: # Double integration
FFTspec = -FFTspec_a*w / (w+EPS)**3
elif times == 1: # Single integration
FFTspec = FFTspec_a*iw / (iw+EPS)**2
elif times == 0: # No integration
FFTspec = FFTspec_a
else:
print 'Error'
FFTspec *= mask # Select frequencies to use
out_w = np.fft.irfft(FFTspec) # Return to time domain
if unwin == True:
out = out_w*window/(window+EPS)**2 # Remove window from time series
else:
out = out_w
if odd == True: # Check for even length time series
return out[:-1] # If not, remove last entry
else:
return out