python中finfo函数_python – 使用scipy fft和ifft以数字方式求解常微分方程

小鱼儿 2023-01-08 08:07 113阅读 0赞

正弦曲线的导数或复指数与其频率成正比,相移π/ 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(Ω).积分几乎是逆操作,但如果存在DC分量,它可能会增加一个冲动：∫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 / 2的采样分量是不可靠的.具体地,每循环仅有2个样本,fs / 2分量的幅度与序列中第一个样本的符号交替.因此,对于实值信号,fs / 2 DFT分量是实值的,具有0或π弧度的相位,即a \* cos(π\* fs \* t \{0,π\}).给定后者,fs / 2分量的导数是-a \*π\* fs \* cos(π\* fs \* t \{π/ 2,3 \*π/ 2\}),对于采样时间t = 0为0 N / FS.

(关于后者,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(Ω)

这是一个单极点低通滤波器.您派生的解决方案有2个问题.

>您没有按fs缩放频率变量w.

>您使用fftfreq,它使用范围-0.5到0.5.你需要-π到π.

实际上你只需要0到π,因为i(t)是实值的.在这

你可以使用rfft和irfft来实现实值信号

跳过计算负频率.

所有这些,你可能仍然对结果感到失望,因为DFT使用信号的周期性扩展.

例子

首先,这是一个简单的例子,1个正弦曲线(以红色绘制),以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讨论了非周期函数的Chebyshev逼近方法(在他的密歇根大学页面免费在线).

如果您在Signal Processing或Mathematics Stack Exchange上询问,您可能会获得更多帮助.