正弦曲线或复指数的导数与其频率成正比，相移π/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)*δ(Ω)
要使用采样序列逼近导数，您必须按采样率 fs 缩放离散时间角频率 ω（范围从 0 到 2π，或 -π 到 π）。例如，假设离散时间频率为 π/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/2 DFT 分量是实值，相位为 0 或 π 弧度，即a*cos(π*fs*t + {0, π})。给定后者，fs/2 分量的导数是-a*π*fs*cos(π*fs*t + {π/2, 3*π/2})，对于采样时间t == n/fs，它是 0。
（关于后者，DFT 的标准三角插值使用余弦，在这种情况下，样本点处的导数将为零。如果同时对信号及其导数进行采样，则不一定正确。采样丢失fs/2 分量相对于信号的相位，但不相对于它的导数。根据您开始采样的时间，fs/2 分量及其导数在采样点处可能不为零。如果通过幸运的是，其中一个在采样时间为 0，另一个将处于峰值，因为它们之间的弧度为 π/2。）
鉴于fs/2 DFT 分量对于实值信号始终为实值，当您在计算导数或积分的过程中将其乘以j 时，这会在结果中引入一个虚值分量.有一个简单的解决方法。如果N 是偶数，只需将索引N/2 处的fs/2 组件清零。另一个问题是除以jΩ 进行积分时除以零。这可以通过在Ω 向量的索引 0 上添加一个小值来解决（例如，finfo(float64).tiny 是最小的双精度浮点数）。
给定Ω = fs * ω，问题中显示的系统在频域中具有以下形式：
H(Ω) = 1 / (g + j*Ω*C)
U(Ω) = I(Ω) * H(Ω)
这是一个单极低通滤波器。您得出的解决方案有两个问题。
- 您没有按
fs 缩放频率变量w。
- 您使用
fftfreq，它使用范围-0.5 到0.5。你需要 -π 到 π。
实际上你只需要 0 到 π 因为i(t) 是实值的。在这个
在这种情况下，您可以将rfft 和irfft 用于实值信号，即
跳过计算负频率。
尽管如此，您可能仍然对结果感到失望，因为 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()
这是一个简单的例子——带宽有限的周期性信号。只需确保以足够高的速率对整数个周期进行采样以避免混叠。
下一个例子是三角波，斜率为 1 和 -1，中心和边缘的导数不连续。理想情况下，导数应该是方波，但完美的计算需要无限带宽。取而代之的是Gibbs ringing 在不连续点周围：
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()
如果您正在求解非周期系统，则 DFT 的隐式周期扩展是有问题的。以下是使用 odeint 和 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()
du/dt 图表覆盖了由diff 估计的导数（红色）与微分方程的计算值（蓝色）。它们在两种情况下大致相等。我将odeint 的初始值设置为u_dft[0]，以表明它返回相同初始值的DFT 解决方案。不同之处在于odeint 解将继续衰减为零，但 DFT 解是周期性的，周期为 2s。在这种情况下，如果对更多 i(t) 进行采样，则 DFT 结果看起来会更好，因为 i(t) 从零开始并衰减到零。
零填充与 DFT 一起用于执行线性卷积。特别是在这种情况下，输入的零填充将有助于将零状态响应的瞬态与其稳态分开。然而，更常用的是拉普拉斯或 z 变换系统函数来分析 ZSR/ZIR。 scipy.signal 拥有分析 LTI 系统的工具，包括以多项式、零极点增益和状态空间形式将连续时间映射到离散时间。
John P. Boyd 在Chebyshev and Fourier Spectral Methods 中讨论了非周期函数的切比雪夫近似方法（在他的密歇根大学页面上免费在线）。
如果您在 Signal Processing 或 Mathematics Stack Exchanges 上提问，您可能会在此类问题上获得更多帮助。