【鼎革‧革鼎】︰ Raspbian Stretch 《六之 J.3‧MIR-10上 》

縱然已知向量

\vec{s_k}, \ k = 0 \cdots N-1

= \left( 1, e^{j \frac{2 \pi \cdot k}{N} \cdot 1}, e^{j \frac{2 \pi \cdot k}{N} \cdot 2}, \cdots , e^{j \frac{2 \pi \cdot k}{N} \cdot (N-1)} \right)

構成了 N離散傅立葉變換 DFT 之正交基底。故該空間裡任何向量都可藉此表現,無有疑焉?想那 N 個角頻率 2 \pi \cdot \frac{k}{N} \cdot f_s 雖多,並不能含括所有頻率乎??比方說這一 2 \pi \cdot \frac{\sqrt{2}}{2} \cdot f_s 不在其中也!

那麼我們將如何詮釋 DFT 的計算結果耶!!

 

※ 參考

scipy.fftpack.fft

scipy.fftpack.fft(x, n=None, axis=-1, overwrite_x=False)
Return discrete Fourier transform of real or complex sequence.

The returned complex array contains y(0), y(1),..., y(n-1) where

y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum().

Parameters:

x : array_like

Array to Fourier transform.

n : int, optional

Length of the Fourier transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].

axis : int, optional

Axis along which the fft’s are computed; the default is over the last axis (i.e., axis=-1).

overwrite_x : bool, optional

If True, the contents of x can be destroyed; the default is False.

Returns:

z : complex ndarray

with the elements:

[y(0),y(1),..,y(n/2),y(1-n/2),...,y(-1)]        if n is even
[y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)]  if n is odd

where:

y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n), j = 0..n-1

See also

ifft
Inverse FFT
rfft
FFT of a real sequence

Notes

The packing of the result is “standard”: If A = fft(a, n), then A[0] contains the zero-frequency term, A[1:n/2] contains the positive-frequency terms, and A[n/2:] contains the negative-frequency terms, in order of decreasingly negative frequency. So for an 8-point transform, the frequencies of the result are [0, 1, 2, 3, -4, -3, -2, -1]. To rearrange the fft output so that the zero-frequency component is centered, like [-4, -3, -2, -1, 0, 1, 2, 3], use fftshift.

Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.

This function is most efficient when n is a power of two, and least efficient when n is prime.

Note that if x is real-valued then A[j] == A[n-j].conjugate(). If x is real-valued and n is even then A[n/2] is real.

If the data type of x is real, a “real FFT” algorithm is automatically used, which roughly halves the computation time. To increase efficiency a little further, use rfft, which does the same calculation, but only outputs half of the symmetrical spectrum. If the data is both real and symmetrical, the dct can again double the efficiency, by generating half of the spectrum from half of the signal.

 

因此 JULIUS O. SMITH III 先生特論

Frequencies in the “Cracks” 》 以解惑勒。

 

正好趁機對照振動現象物理說明哩◎

依據牛頓第二運動定律,一個受驅振子的方程式為

F(t)-kx-c\frac{\mathrm{d}x}{\mathrm{d}t}=m\frac{\mathrm{d}^2x}{\mathrm{d}t^2},一般將之改寫為

\frac{\mathrm{d}^2x}{\mathrm{d}t^2} + 2\zeta\omega_0\frac{\mathrm{d}x}{\mathrm{d}t} + \omega_0^2 x = \frac{F(t)}{m},此處 \omega_0 = \sqrt{\frac{k}{m}}  稱為『無阻尼』角頻率,而 \zeta = \frac{c}{2 \sqrt{mk}} 叫做『阻尼比率』。如果『外力F(t) = 0,那個方程式就成了『阻尼振子』的方程式

\frac{\mathrm{d}^2x}{\mathrm{d}t^2} + 2\zeta\omega_0\frac{\mathrm{d}x}{\mathrm{d}t} + \omega_0^2 x = 0,當 \zeta \leq 1 它的解是

x(t) = A \mathrm{e}^{-\zeta \omega_0 t} \ \sin \left( \sqrt{1-\zeta^2} \ \omega_0 t + \phi \right),此處 \phi 是相位角。

假使 \zeta > 1 它的解是

x(t) = C_1 \ \mathrm{e}^{\left(- \zeta + \sqrt{\zeta^2 -1}\right) \omega_0 t} + C_2  \ \mathrm{e}^{\left(- \zeta - \sqrt{\zeta^2 -1}\right) \omega_0 t}

如果我們從 \zeta 的值來看『阻尼振子』的系統行為,當 \zeta > 1 時,這一個系統已經『振動』不起來了,通常叫做『過阻尼』,負數的『指數項』使得系統的能量隨時間逐漸減少,\zeta 的數值愈大能量減少將慢愈遲回到平衡。當 \zeta = 1 時,這一個系統也『振動』不起來了,通常稱之為『臨界阻尼』,此時系統會用最快的方式設法回到平衡,這個可是『關門』系統的『最佳解』!!。當 \zeta < 1 時,這樣的諧振子系統稱作『低阻尼』,這時系統用著『低於無阻尼』的『頻率』振動,系統的『振幅』隨著時間以 \mathrm{e}^{-\zeta \omega_0 t} 為比率逐漸減小以至於『不振動 』為止。事實上從自然界中來的一般現象都會比『理論值』更快的到達『停止點』,就像說不只有『動摩擦力』與『靜摩擦力』之區分,摩擦力的『速度相關性』也不是這麼『簡單的正比』之假設,更別說理論上還有著『摩擦生熱』的問題必須要考慮。我們也許可以說為著追求『基本現象的理解』,通常都會『假設』了一些數學上『解答問題』的『理想條件』。

現在談談受外力影響下的受驅振子︰

階躍 Step  外力

假設此系統的 \zeta < 1,初始位置 x_0 = 0,在 t = 0^+ 時受到如下的階躍外力︰

{F(t) \over m} = \begin{cases} \omega _0^2 & t \geq 0 \\ 0 & t < 0 \end{cases}

Dirac_distribution_CDF.svg

它的解是

x(t) = 1 - \mathrm{e}^{-\zeta \omega_0 t} \frac{\sin \left( \sqrt{1-\zeta^2} \ \omega_0 t + \varphi \right)}{\sin(\varphi)},此處相位角 \varphi\cos \varphi = \zeta 所決定。

這個系統因為零點時刻突然受到固定大小的外力 m  \ {\omega_0}^2 所驅動,震盪以 \mathrm{e}^{-\zeta \omega_0 t} 為比率逐漸增大,一般用 \tau = \frac{1}{\zeta \omega_0} 為時間尺度來衡量這個變化,每一 \tau 單位時間,系統將以 \mathrm{e}^{-1} 為比率改變振幅,在物理上稱之為『弛豫時間』Relaxation Time,工程上常用多的 \tau 單位時間,來談震盪達到預期大小的『安定時間』settling time。果真是『風吹枝擺 』,待其風歇『搖曳而止』!!

頻率為 \omega 的正弦驅動力

此時系統的方程式為

\frac{\mathrm{d}^2x}{\mathrm{d}t^2} + 2\zeta\omega_0\frac{\mathrm{d}x}{\mathrm{d}t} + \omega_0^2 x = \frac{1}{m} F_0 \sin(\omega t)

220px-Sin.svg

F_0 是驅動力的振幅大小。在線性微分方程式如 \hat{L} x(t) = F(t) 的『求解』裡,如過『\Box』是 \hat{L} x(t) = 0 的一個解,『\bigcirc』是 \hat{L} x(t) = F(t) 一個『特解』,那麼『c \ \Box +  \ \bigcirc』就是該方程式的『通解』。我們已經知道 F(t) = 0 的『低阻尼振子』之解在若干個弛豫時間後數值將變得太小了,所以它對於系統長時間之後的『行為 』沒有太多的貢獻。因此我們說這個系統的『穩態解』steady-state solution 是

x(t) = \frac{F_0}{m Z_m \omega} \sin(\omega t + \phi),此處

Z_m = \sqrt{\left(2\omega_0\zeta\right)^2 + \frac{1}{\omega^2}\left(\omega_0^2 - \omega^2\right)^2}

是『響應阻抗』函數。而 \phi 是驅動力引發的相位角,可由

\phi = \arctan\left(\frac{\omega_0^2-\omega^2}{2\omega \omega_0\zeta}\right)

所決定,一般它表達著相位『遲滯』 lag 現象。

300px-Resonance

如果使用不同頻率 \omega 的驅動力,當\omega = \omega_r = \omega_0\sqrt{1-2\zeta^2} 時,系統的響應振幅最大,這稱之為『共振』resonant,這一個頻率就叫做『共振頻率』。

請參考左圖。

物理上所說的『慣性』是指一個系統遭受外力時,它會發生『抵抗變化』的作為。這就是『響應阻抗』和『相位遲滯』的物理原由與命名由來。假使考察穩態解,我們是否可以講︰『原因 』── F_0 \sin(\omega t) ── 產生成正比之『結果』── x(t) = \frac{F_0}{m Z_m \omega} \sin(\omega t + \phi) ── 的呢??

─── 摘自《【Sonic π】聲波之傳播原理︰振動篇