6

I have an ordinary differential equation in time domain as follows: 

C*du/dt = -g*u + I

where I = A*t/tau*exp^(1-t/tau)

in the freq domain:

u(w) = I(w)/(g*(1+C/g*j*w))

j being the complex number sqrt(-1)

hence i can get u(t) by going into the freq domain using fast Fourier transform (fft) and then back using ifft.

the codes:

t = np.linspace(0.,499.9,5000)
I = q*t*np.exp(1-t/Tau_ca)/Tau_ca
u1 = np.fft.ifft(np.fft.fft(I)/(G_l*(1.+1.j*(C_m/G_l)*np.fft.fftfreq(t.shape[-1]))))

However, when I compare this u(t) obtained with other methods such as numerical integration of the differential equation or its analytical form, it is not correct. I have tried and been unsuccessful at figuring out where my mistakes are.

please enlighten.

4

1 回答 1

11

正弦曲线或复指数的导数与其频率成正比，相移π/2. 对于复指数，相移相当于乘以j。例如，d/dt exp(j*Ω*t) == j*Ω * exp(j*Ω*t) == Ω * exp(j*π/2) * exp(j*Ω*t) == Ω * exp(j*(Ω*t + π/2))。所以如果你有傅里叶变换对u(t) <=> U(Ω)，那么du/dt <=> jΩ * U(Ω)。积分几乎是逆运算，但如果有直流分量，它可能会增加脉冲：∫udt <=> U(Ω) / jΩ + π*U(0)*δ(Ω)

要使用采样序列逼近导数，您必须通过采样率来缩放离散时间角频率 ω（范围从 0 到 2π，或 -π 到 π）fs。例如，假设离散时间频率为 π/2 弧度/样本，例如序列[1, 0, -1, 0, ...]。在原始信号中，这对应于fs/4。导数是d/dt cos(2*π*fs/4*t) == d/dt cos(π/2*fs*t) == -π/2*fs*sin(π/2*fs*t) == π/2*fs*cos(π/2*fs*t + π/2)。

您必须以fs信号带宽的两倍以上进行采样。准确的采样分量 fs/2是不可靠的。具体来说，每个周期只有 2 个样本，fs/2分量的幅度会交替序列中第一个样本的符号。因此，对于实值信号，fs/2DFT 分量是实值，相位为 0 或 π 弧度，即a*cos(π*fs*t + {0, π})。给定后者，fs/2分量的导数为-a*π*fs*cos(π*fs*t + {π/2, 3*π/2})，对于采样时间为 0 t == n/fs。

（关于后者，DFT 的标准三角插值使用余弦，在这种情况下，样本点处的导数将为零。如果同时对信号及其导数进行采样，则不一定正确。采样会丢失相位分量相对于fs/2信号，但不相对于它的导数。根据您开始采样的时间，fs/2分量及其导数在采样点可能都不为零。如果幸运的话，其中一个在采样点为 0有时，另一个会达到顶峰，因为它们之间的π/2弧度相差很大。）

鉴于fs/2DFT 分量对于实值信号始终是实值，当您在计算导数或积分的过程中将其乘以时j，这会在结果中引入虚值分量。有一个简单的解决方法。如果是偶数，只需将index 处N的组件清零。另一个问题是除以进行积分时除以零。这可以通过向向量的索引 0 添加一个小值来解决（例如，最小的双精度浮点数）。fs/2N/2jΩΩfinfo(float64).tiny

给定Ω = fs * ω，问题中显示的系统在频域中具有以下形式：

H(Ω) = 1 / (g + j*Ω*C)
U(Ω) = I(Ω) * H(Ω)

它是一个单极低通滤波器。您得出的解决方案有两个问题。

  1. 您没有按比例缩放频率w变量fs。
  2. 您使用fftfreq，它使用范围 -0.5 到 0.5。你需要 -π 到 π。实际上你只需要 0 到 π，因为i(t)它是实值的。在这种情况下，您可以将rfftandirfft用于实值信号，从而跳过计算负频率。

尽管如此，您可能仍然对结果感到失望，因为 DFT 使用信号的周期性扩展。

例子

首先，这是一个 1 Hz 正弦曲线的简单示例（以红色绘制），以 1024 个样本/秒的速度采样 2 秒，并通过 DFT 计算其导数（以蓝色绘制）：

from pylab import *

fs = 1024
t = arange(2*fs, dtype=float64) / fs
N = len(t) // 2 + 1    # non-negative frequencies
w = 2 * pi / len(t) * arange(N)
Omega = w * fs

x0 = cos(2*pi*t)    # w0 == 2*pi/fs
X0 = rfft(x0);
# Zero the fs/2 component. It's zero here anyway.
X0[-1] = 0  
dx0 = irfft(X0 * 1j*Omega)
plot(t, x0, 'r', t, dx0, 'b')
show()

dt0

这是一个简单的例子——带宽有限的周期性信号。只需确保以足够高的速率对整数个周期进行采样以避免混叠。

下一个例子是三角波，斜率为 1 和 -1，中心和边缘的导数不连续。理想情况下，导数应该是方波，但完美的计算需要无限带宽。相反，吉布斯在不连续性周围响起：

t2 = t[:fs]
m = len(t) // (2*len(t2))
x1 = hstack([t2, 1.0 - t2] * m)
X1 = rfft(x1)
X1[-1] = 0
dx1 = irfft(X1 * 1j*Omega)
plot(t, x1, 'r', t, dx1, 'b')
show()

dt1

如果您正在求解非周期性系统，则 DFT 的隐式周期性扩展是有问题的。这是使用odeintDFT 和 DFT（tau设置为 0.5s；g并C设置为 1 Hz 转角频率）的问题系统的解决方案：

from scipy.integrate import odeint

a = 1.0; tau = 0.5
g = 1.0; C = 1.0 / (2 * pi)
i = lambda t: a * t / tau * exp(1 - t / tau)
f = lambda u, t: (-g*u + i(t)) / C

H = 1 / (g + 1j*Omega*C)  # system function
I = rfft(i(t))
I[-1] = 0
U_DFT = I * H
u_dft = irfft(U_DFT)
u_ode = odeint(f, u_dft[0], t)[:,0]

td = t[:-1] + 0.5/fs
subplot('221'); title('odeint u(t)');
plot(t, u_ode)
subplot('222'); title('DFT u(t)');
plot(t, u_dft)
subplot('223'); title('odeint du/dt')
plot(td, diff(u_ode)*fs, 'r',
     t, (-g*u_ode + i(t)) / C, 'b')           
subplot('224'); title('DFT du/dt')
plot(td, diff(u_dft)*fs, 'r',
     t, (-g*u_dft + i(t)) / C, 'b')
show()

dt2

du/dt图表覆盖了由（红色）估计的导数与diff微分方程（蓝色）的计算值。它们在两种情况下大致相等。我将初始值设置为odeinttou_dft[0]以表明它返回相同初始值的 DFT 解决方案。不同之处在于odeint解会继续衰减到零，但 DFT 解是周期性的，周期为 2s。在这种情况下，如果对更多 i(t) 进行采样，则 DFT 结果看起来会更好，因为 i(t) 从零开始并衰减到零。

DFT 使用零填充来执行线性卷积。特别是在这种情况下，输入的零填充将有助于将零状态响应的瞬态与其稳态分开。然而，更常用的是拉普拉斯或 z 变换系统函数来分析 ZSR/ZIR。scipy.signal拥有分析 LTI 系统的工具，包括以多项式、零极点增益和状态空间形式将连续时间映射到离散时间。

John P. Boyd 在Chebyshev 和 Fourier Spectral Methods中讨论了非周期函数的 Chebyshev 逼近方法（在他的密歇根大学页面上免费在线）。

如果您在信号处理或数学堆栈交换中提问，您可能会在此类问题上获得更多帮助。

于 2012-11-29T22:23:07.253 回答