Considerations in the design of a ShortTimeFFT window?
I was very excited by the new ShorftTimeFFT class, and am beginning to introduce it into my workflow. My team is interested in audio simulation and source waveform reconstruction from microphone arrays, and the inverse STFT functionality is a big draw to the new class. Our first use case is an integrator and differentiator. We implement this with a forward/inverse stft pair, and the frequency domain calculus definitions. We were surprised that the choice of window changes the derivative result in particular, at a -50 dB level. Unexpectedly, it seems that lower order windows work much better than higher order. For example, the hann window computes a better inverse stft than the nuttall window. We would like to get a better sense of why and get better insight into this new class. The ShortTImeFFT example uses a Gaussian window, which is an unusual window choice in my experience. Is there a good reference that can guide the selection of a window with the application of an inverse stft in mind? Thank you, Ned Richards
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
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.Brunn@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
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$")
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:
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
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 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.Brunn@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-tran...
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.Brunn@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
participants (2)
-
Dietrich Brunn -
Edward Richards