The example got truncated while pasting - here is the remainder:
...
ax1.plot(t, dx, alpha=0.5, label="STFT-based $dx/dt$")
for ax_ in (*axx0, ax1):
ax_.legend()
ax_.grid(True)
plt.show()
Cheers, dietrich
On Monday, 4 March 2024 16:51:51 CET Dietrich Brunn wrote:
> 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:
>
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
