Hi Ned, I took the time to think about this a little bit: Unfortunately I am not aware of useful literature for choosing STFT windows. Note that the concept of a dual window for the ISTFT does not seem to be widely known in the signal processing community. The only other implementation I know, which provides a dual window is LTFAT [1]
Designing a STFT based differentiator is a bit involving, so let's discuss an
FFT based one
first: The continuous-time signal representation of a sampled signal of finite
length can be
expressed by a complex-valued Fourier series E. g., for signal of duration
one, we can
write
x(t) = Sum[ X[l] exp(2jπ(l Δ f) t ) ] with l being the summation index.
The coefficients X[l] can be calculated with an FFT ([2] discusses this from a
different
angle). Differentiation with respect to time gives
d/dt x(t) = Sum[ X[l] 2jπ(l Δ f) exp(2jπ(l Δ f) t) ] .
Note that x(t) is assumed to be periodic. Hence a discontinuity between start
and end of
signal produces ringing due to Gibb's phenomenon.
The key insight of the ShortTimeFFT implementation is that any signal can be
represented
by a series expansion of time- and frequency-shifted dual windows, i.e.,
x(t) = Sum[ S[q, p] d(t - p Δt) exp(2jπ(q Δ f) t ) ]
with (p Δt) representing the time shift, (q Δ f) the frequency shift, d(t) the
dual window and
S[p, q] the STFT coefficient (consult [3] for details).
Differentiation with respect to time gives
d/dt x(t) = Sum[ S[q, p] 2jπ(q Δ f) d(t - p Δt) exp(2jπ(q Δ f) t ) ] +
S[q, p] exp(2jπ(q Δ f) t ) d/dt d(t - p Δt) ]
Note that:
* The window dependent term is parameterized by the derivative of the dual
window d/dt
d(t - p Δt).
* If the signal is not periodic in each slice which is stenciled out by the
sliding window,
Gibb's phenomenon will strike again (this is what you probably observe in your
simulation). Hence it is a good idea to choose a dual window, like the Hann
window, which
suppresses discontinuities at the beginning and end of the signal slice.
* Not only the dual window but also its derivative needs to suppress those
discontinuities.
The following example shows an example of STFT based differentiation. The hop
width is
chosen small enough to supress Gibb's phenomenon in the dual window derivative.
Note
that Gibb's phenomenon can be observed at the beginning and the end of the
signal.
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import ShortTimeFFT
from scipy.signal import windows
from scipy.fft import rfft, rfftfreq, irfft
# Create periodic test signal and its derivative:
n, T = 1000, 1/1000 # samples and sampling interval for 1 second signal
t = np.arange(n) * T # time stamps
# Create single frequency signal and derivative:
f = rfftfreq(n, T)
k_c = f.searchsorted(7)
omega_c = 2*np.pi*f[k_c]
X = np.zeros(len(f))
X[k_c] = n/2/omega_c
x = irfft(X, n=n)
y = irfft(2j*np.pi*f*X, n=n) # dx / dt
# Two ShortTimeFFT instances are needed:
kw = dict(hop=10, fs=1/T,fft_mode='onesided', phase_shift=None)
SFT = ShortTimeFFT.from_dual(windows.hann(32, sym=False), **kw)
# Differentiate dual window per FFT:
diff_dual_win = irfft(rfft(SFT.dual_win) * 2j*np.pi*rfftfreq(SFT.m_num))
dSFT = ShortTimeFFT.from_dual(diff_dual_win, **kw)
# Perform the filtering:
S_x = SFT.stft(x)
dS_x = 2j*np.pi*SFT.f[:, np.newaxis] * S_x
dx = SFT.istft(dS_x, k1=n) + dSFT.istft(S_x, k1=n)
# Plot windows:
fg0, axx0 = plt.subplots(2, 1, sharex='all', tight_layout=True)
axx0[0].set(title="Windows")
axx0[1].set(title="Dual Windows")
axx0[0].plot(SFT.win, '.-', alpha=0.5, label='SFT')
axx0[0].plot(dSFT.win, '.-', alpha=0.5, label='dSFT')
axx0[1].plot(SFT.dual_win, '.-', alpha=0.5, label='SFT')
axx0[1].plot(dSFT.dual_win, '.-', alpha=0.5, label='dSFT')
# Plot signal:
fg1, ax1 = plt.subplots()
ax1.set(title=rf"STFT-based Differentiator for $f_c={f[k_c]}\,$Hz Signal",
xlabel="Time $t$", ylabel="Amplitude")
ax1.plot(t, omega_c * x, alpha=0.5, label=r"$x(t) \cdot 2\pi f_c$")
ax1.plot(t, y, '--', alpha=0.5, label="$y = dx/dt$")
signature.asc
Description: This is a digitally signed message part.
_______________________________________________ SciPy-Dev mailing list -- scipy-dev@python.org To unsubscribe send an email to scipy-dev-le...@python.org https://mail.python.org/mailman3/lists/scipy-dev.python.org/ Member address: arch...@mail-archive.com
