Dietrich:
Thank you for taking the time to both derive a result and make an example.
I am going to keep mulling this over for a bit, but I wanted to say thank
you for your help.
I found the expression
> x(t) = Sum[ S[q, p] d(t - p Δt) exp(2jπ(q Δ f) t ) ]
very useful, and I appreciate your highlighting it for me. The notion that
the dual window is sampling across the chunks helped me see that making the
hop size smaller could improve the inverse stft result.
I was using the COLA spacing of window length / 4 for the nuttall window,
and when I switched to the hann window this same hop length was 2x higher
resolution than the first COLA hop length. Finally, I changed the nuttall
hop length to window length / 8, which gave very good reconstruction
results (away from the edges).
Best,
Ned
On Mon, Mar 4, 2024 at 7:51 AM Dietrich Brunn <dietrich.br...@web.de> 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:
>
> 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$")
> 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
>
>
>
> [1] https://ltfat.org/doc/gabor/
>
> [2] https://scipy.github.io/devdocs/tutorial/signal.html#spectral-analysis
>
> [3]
> https://scipy.github.io/devdocs/tutorial/signal.html#short-time-fourier-transform
>
> On Saturday, 2 March 2024 22:10:45 CET Edward Richards wrote:
>
> > Thank you for your response, Dietrich.
>
> >
>
> > We are currently implementing the differentiation as you suggested: (1)
>
> > STFT the signal, (2) multiply each slice by j2πf, (3) ISTFT. We are not
>
> > using a rectangular window, though I can see how it performs. I default
> to
>
> > a nuttall for DSP, and its relatively poor performance lead to my
> question.
>
> > Is there a way of predicting if an operation does depend on the window?
> If
>
> > so, is there a standard practice for selection?
>
> >
>
> > This procedure does show promise. The integration step is much simpler
> with
>
> > the STFT than a simple FFT/IFFT.
>
> >
>
> > On Sat, Mar 2, 2024 at 1:49 AM Dietrich Brunn <dietrich.br...@web.de>
> wrote:
>
> > > Hi Ned,
>
> > >
>
> > > it is nice to see that the new STFT functionality is actually being
> used.
>
> > >
>
> > > > Our first use case is an integrator and differentiator. We implement
>
> > > > this
>
> > > >
>
> > > > with a forward/inverse stft pair, and the frequency domain calculus
>
> > > >
>
> > > > definitions.
>
> > >
>
> > > May I ask how you implemented the differentiation in the STFT space?
> Note
>
> > > that differentiation *is* window dependent. An ad-hoc approach would be
>
> > > using a rectangular dual window (utilizing `ShortTimeFFT.from_dual`)
> and
>
> > > multiplying the STFT by j2πf.
>
> > >
>
> > > Cheers, dietrich
>
>
>
_______________________________________________
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
