6

次のようなordinary differential equationインタイムドメインがあります。

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

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

周波数ドメイン:

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

j複素数であることsqrt(-1)

したがって、高速フーリエ変換( )u(t)を使用して周波数領域に入り、 を使用し て戻ることができます。fftifft

コード:

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]))))

しかし、これu(t)を微分方程式の数値積分やその解析形式などの他の方法と比較すると、正しくありません。私は自分の間違いがどこにあるかを理解しようとしましたが、うまくいきませんでした。

啓発してください。

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(Ω). 積分はほぼ逆の操作ですが、DC 成分がある場合はインパルスが追加される可能性があります。∫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信号の帯域幅の 2 倍以上でサンプリングする必要があります。正確にサンプリングされたコンポーネント fs/2は信頼できません。具体的には、1 サイクルあたり 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成分とその導関数の両方がサンプル点で非ゼロになる場合があります. 運が良ければ、そのうちの 1 つがサンプルで 0 の場合π/2ラジアンが離れているため、もう一方はピークになります。)

fs/2実数値の信号に対して DFT コンポーネントが常に実数値であることを考えるとj、導関数または積分を計算する過程で DFT コンポーネントを乗算すると、結果に虚数値のコンポーネントが導入されます。簡単な回避策があります。が偶数の場合は、 indexNのコンポーネントをゼロにします。もう 1 つの問題は、積分のために で割るときの 0 による除算です。これは、ベクトルのインデックス 0 に小さな値を追加することで解決できます(たとえば、最小の倍精度浮動小数点数です)。fs/2N/2jΩΩfinfo(float64).tiny

が与えられΩ = fs * ωた場合、問題に示されているシステムは、周波数領域で次の形式になります。

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

単極ローパスフィルターです。導き出した解には 2 つの問題があります。

  1. w周波数変数を でスケーリングしていませんfs。
  2. fftfreq-0.5 から 0.5 の範囲を使用する を使用します。-π から π が必要です。i(t)は実数値なので、実際には 0 から π だけが必要です。この場合、負の周波数の計算をスキップする実数値信号にrfftandを使用できます。irfft

とはいえ、DFT は信号の周期的な拡張を使用するため、結果にまだ失望する可能性があります。

例

まず、1024 サンプル/秒で 2 秒間サンプリングされた 1 Hz の正弦波 (赤でプロット) と、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 の暗黙の周期的拡張は問題になります。odeintと DFT の両方を使用した問題のシステムのソリューションを次に示します(tauは 0.5 秒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。どちらの場合もほぼ同じです。同じ初期値に対して DFT 解を返すことを示すために、 odeinttoの初期値を設定します。u_dft[0]違いは、odeint解がゼロまで減衰し続けることですが、DFT 解は 2 秒周期で周期的です。この場合、i(t) がゼロから始まりゼロに減衰するため、より多くの i(t) がサンプリングされると、DFT の結果はより良く見えます。

DFT でゼロ パディングを使用して線形畳み込みを実行します。特にこの場合、入力のゼロ パディングは、ゼロ状態応答の過渡状態をその定常状態から分離するのに役立ちます。ただし、より一般的には、ZSR/ZIR を分析するためにラプラスまたは z 変換システム関数が使用されます。scipy.signalには、連続時間から離散時間への多項式、零点-極-ゲイン、状態空間形式へのマッピングなど、LTI システムを解析するためのツールがあります。

John P. Boyd が、 Chebyshev および Fourier Spectral Methodsで非周期関数の Chebyshev 近似法について説明しています (ミシガン大学の彼のページでオンラインで無料)。

このような質問については、 Signal ProcessingまたはMathematics Stack Exchangesで質問すると、おそらくより多くの助けが得られるでしょう。

于 2012-11-29T22:23:07.253 に答える