Re: Considerations in the design of a ShortTimeFFT window?
Edward Richards <[email protected]> Wed, 6 Mar 2024 10:52:34 -0800
| Newsgroups | gmane.comp.python.scientific.devel |
|---|---|
| Message-ID | <CAADEu5YRT-hDGk7F4NCmmNDgYftM8Umo5K0Gedqq31MOpeYCUQ@mail.gmail.com> |
--===============1865187239829479041== Content-Type: multipart/alternative; boundary="000000000000d8cac40613027684" --000000000000d8cac40613027684 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable 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) =3D Sum[ S[q, p] d(t - p =CE=94t) exp(2j=CF=80(q =CE=94 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=E2=80=AFAM Dietrich Brunn <[email protected]= e> 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 conce= pt > 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 discus= s > 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) =3D Sum[ X[l] exp(2j=CF=80(l =CE=94 f) t ) ] with l being the summ= ation 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) =3D Sum[ X[l] 2j=CF=80(l =CE=94 f) exp(2j=CF=80(l =CE=94 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) =3D Sum[ S[q, p] d(t - p =CE=94t) exp(2j=CF=80(q =CE=94 f) t ) ] > > with (p =CE=94t) representing the time shift, (q =CE=94 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) =3D Sum[ S[q, p] 2j=CF=80(q =CE=94 f) d(t - p =CE=94t) exp(2j= =CF=80(q =CE=94 f) t ) ] + > > S[q, p] exp(2j=CF=80(q =CE=94 f) t ) d/dt d(= t - p =CE=94t) ] > > Note that: > > * The window dependent term is parameterized by the derivative of the > dual window d/dt d(t - p =CE=94t). > > * 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 th= e > beginning and end of the signal slice. > > * Not only the dual window but also its derivative needs to suppress thos= e > 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 =3D 1000, 1/1000 # samples and sampling interval for 1 second sign= al > > t =3D np.arange(n) * T # time stamps > # Create single frequency signal and derivative: > > f =3D rfftfreq(n, T) > k_c =3D f.searchsorted(7) > omega_c =3D 2*np.pi*f[k_c] > X =3D np.zeros(len(f)) > X[k_c] =3D n/2/omega_c > x =3D irfft(X, n=3Dn) > y =3D irfft(2j*np.pi*f*X, n=3Dn) # dx / dt > > > > # Two ShortTimeFFT instances are needed: > > kw =3D dict(hop=3D10, fs=3D1/T,fft_mode=3D'onesided', phase_shift=3DNone= ) > SFT =3D ShortTimeFFT.from_dual(windows.hann(32, sym=3DFalse), **kw) > > # Differentiate dual window per FFT: > > diff_dual_win =3D irfft(rfft(SFT.dual_win) * 2j*np.pi*rfftfreq(SFT.m_num= )) > dSFT =3D ShortTimeFFT.from_dual(diff_dual_win, **kw) > > > > # Perform the filtering: > > S_x =3D SFT.stft(x) > dS_x =3D 2j*np.pi*SFT.f[:, np.newaxis] * S_x > dx =3D SFT.istft(dS_x, k1=3Dn) + dSFT.istft(S_x, k1=3Dn) > > > > # Plot windows: > > fg0, axx0 =3D plt.subplots(2, 1, sharex=3D'all', tight_layout=3DTrue) > axx0[0].set(title=3D"Windows") > axx0[1].set(title=3D"Dual Windows") > axx0[0].plot(SFT.win, '.-', alpha=3D0.5, label=3D'SFT') > axx0[0].plot(dSFT.win, '.-', alpha=3D0.5, label=3D'dSFT') > axx0[1].plot(SFT.dual_win, '.-', alpha=3D0.5, label=3D'SFT') > axx0[1].plot(dSFT.dual_win, '.-', alpha=3D0.5, label=3D'dSFT') > > > > # Plot signal: > > fg1, ax1 =3D plt.subplots() > ax1.set(title=3Drf"STFT-based Differentiator for $f_c=3D{f[k_c]}\,$Hz Sig= nal", > xlabel=3D"Time $t$", ylabel=3D"Amplitude") > ax1.plot(t, omega_c * x, alpha=3D0.5, label=3Dr"$x(t) \cdot 2\pi f_c$") > ax1.plot(t, y, '--', alpha=3D0.5, label=3D"$y =3D dx/dt$") > ax1.plot(t, dx, alpha=3D0.5, label=3D"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-analysi= s > > [3] > https://scipy.github.io/devdocs/tutorial/signal.html#short-time-fourier-t= ransform > > 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=CF=80f, (3) ISTFT. We ar= e 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=E2=80=AFAM Dietrich Brunn <Dietrich.Brunn@w= eb.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 implemen= t > > > > > 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=CF=80f. > > > > > > > > Cheers, dietrich > > > --000000000000d8cac40613027684 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr"><div>Dietrich:</div><div>Thank you for taking the time to = both derive a result and make an example. I am going to keep mulling this o= ver for a bit, but I wanted to say thank you for your help.</div><div><br><= /div><div>I found the expression</div><div>> x(t) =3D Sum[=C2=A0 S[q, p]= d(t - p =CE=94t) exp(2j=CF=80(q =CE=94 f) t ) ] <br></div><div>very useful= , and I appreciate your highlighting it for me. The notion that the dual wi= ndow is sampling across the chunks helped me see that making the hop size s= maller could improve the inverse stft result.=C2=A0</div><div><br></div><di= v>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 h= op length to window length / 8, which gave very good reconstruction results= (away from the edges).</div><div><br></div><div>Best,</div><div>Ned<br></d= iv></div><br><div class=3D"gmail_quote"><div dir=3D"ltr" class=3D"gmail_att= r">On Mon, Mar 4, 2024 at 7:51=E2=80=AFAM Dietrich Brunn <<a href=3D"mai= lto:[email protected]">[email protected]</a>> wrote:<br></div><b= lockquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;border-le= ft:1px solid rgb(204,204,204);padding-left:1ex"> <div><p style=3D"margin:0px">Hi Ned,</p> <p style=3D"margin:0px">I took the time to think about this a little bit: U= nfortunately 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 w= idely known in the signal processing community. The only other implementati= on I know, which provides a dual window is LTFAT [1] </p> <br><p style=3D"margin:0px">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 expresse= d by a complex-valued Fourier series=C2=A0 E. g., for signal of duration on= e, we can write</p> <br><p style=3D"margin:0px">x(t) =3D Sum[=C2=A0 X[l]=C2=A0 exp(2j=CF=80(l = =CE=94 f) t ) ] with l being the summation index. </p> <br><p style=3D"margin:0px">The coefficients X[l] can be calculated with an= FFT ([2] discusses this from a different angle). Differentiation with resp= ect to time gives </p> <br><p style=3D"margin:0px">d/dt x(t) =3D Sum[=C2=A0 X[l]=C2=A0 2j=CF=80(l = =CE=94 f) exp(2j=CF=80(l =CE=94 f) t) ] .</p> <br><p style=3D"margin:0px">Note that x(t) is assumed to be periodic. Hence= a discontinuity between start and end of signal produces ringing due to Gi= bb's phenomenon.</p> <br><br><p style=3D"margin:0px">The key insight of the ShortTimeFFT impleme= ntation is that any signal can be represented by a series expansion of time= - and frequency-shifted dual windows, i.e.,</p> <br><p style=3D"margin:0px">x(t) =3D Sum[=C2=A0 S[q, p] d(t - p =CE=94t) ex= p(2j=CF=80(q =CE=94 f) t ) ] </p> <br><p style=3D"margin:0px">with (p =CE=94t) representing the time shift,= =C2=A0 (q =CE=94 f) the frequency shift, d(t) the dual window and S[p, q] t= he STFT coefficient (consult [3] for details).</p> <p style=3D"margin:0px">Differentiation with respect to time gives </p> <br><p style=3D"margin:0px">d/dt x(t) =3D Sum[ S[q, p] 2j=CF=80(q =CE=94 f)= =C2=A0 d(t - p =CE=94t) exp(2j=CF=80(q =CE=94 f) t ) ]=C2=A0 + </p> <p style=3D"margin:0px">=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2= =A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0= =C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 S[q, p] exp(2j=CF=80(q =CE=94 f) t )= =C2=A0 d/dt d(t - p =CE=94t) ] </p> <br><p style=3D"margin:0px">Note that:</p> <br><p style=3D"margin:0px">* The window dependent term is parameterized by= the derivative of the=C2=A0 dual window d/dt d(t - p =CE=94t).</p> <br><p style=3D"margin:0px">* If the signal is not periodic in each slice w= hich is stenciled out by the sliding window, Gibb's phenomenon will str= ike 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 suppres= ses discontinuities at the beginning and end of the signal slice.=C2=A0 </p= > <br><p style=3D"margin:0px">* Not only the dual window but also its derivat= ive needs to suppress those discontinuities. </p> <br><p style=3D"margin:0px">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=C2=A0 that Gibb's= phenomenon=C2=A0 can be observed at the beginning and the end of the signa= l.</p> <br><p style=3D"margin:0px"><span style=3D"font-family:JetBrains Mono"><spa= n style=3D"color:rgb(207,142,109)">import </span><span style=3D"color:rgb(1= 88,190,196)">matplotlib.pyplot </span><span style=3D"color:rgb(207,142,109)= ">as </span><span style=3D"color:rgb(188,190,196)">plt</span></span></p> <p>=C2=A0<span style=3D"color:rgb(207,142,109)">import </span><span style= =3D"color:rgb(188,190,196)">numpy </span><span style=3D"color:rgb(207,142,1= 09)">as </span><span style=3D"color:rgb(188,190,196)">np</span></p> <p>=C2=A0<span style=3D"color:rgb(207,142,109)">from </span><span style=3D"= color:rgb(188,190,196)">scipy.signal </span><span style=3D"color:rgb(207,14= 2,109)">import </span><span style=3D"color:rgb(188,190,196)">ShortTimeFFT</= span></p> <p>=C2=A0<span style=3D"color:rgb(207,142,109)">from </span><span style=3D"= color:rgb(188,190,196)">scipy.signal </span><span style=3D"color:rgb(207,14= 2,109)">import </span><span style=3D"color:rgb(188,190,196)">windows</span>= </p> <p>=C2=A0<span style=3D"color:rgb(207,142,109)">from </span><span style=3D"= color:rgb(188,190,196)">scipy.fft </span><span style=3D"color:rgb(207,142,1= 09)">import </span><span style=3D"color:rgb(188,190,196)">rfft, rfftfreq, i= rfft</span></p> <p>=C2=A0</p><p>=C2=A0<span style=3D"color:rgb(122,126,133)"># Create perio= dic test signal and its derivative:</span></p> <p>=C2=A0<span style=3D"color:rgb(188,190,196)">n, T =3D </span><span style= =3D"color:rgb(42,172,184)">1000</span><span style=3D"color:rgb(188,190,196)= ">, </span><span style=3D"color:rgb(42,172,184)">1</span><span style=3D"col= or:rgb(188,190,196)">/</span><span style=3D"color:rgb(42,172,184)">1000=C2= =A0 </span><span style=3D"color:rgb(122,126,133)"># samples and sampling in= terval for 1 second signal</span></p> <p>=C2=A0<span style=3D"color:rgb(188,190,196)">t =3D np.arange(n) * T=C2= =A0 </span><span style=3D"color:rgb(122,126,133)"># time stamps<br># Create= single frequency signal and derivative:</span></p> <p>=C2=A0<span style=3D"color:rgb(188,190,196)">f =3D rfftfreq(n, T)<br>k_c= =3D f.searchsorted(</span><span style=3D"color:rgb(42,172,184)">7</span><s= pan style=3D"color:rgb(188,190,196)">)<br>omega_c =3D </span><span style=3D= "color:rgb(42,172,184)">2</span><span style=3D"color:rgb(188,190,196)">*np.= pi*f[k_c]<br>X =3D np.zeros(</span><span style=3D"color:rgb(136,136,198)">l= en</span><span style=3D"color:rgb(188,190,196)">(f))<br>X[k_c] =3D n/</span= ><span style=3D"color:rgb(42,172,184)">2</span><span style=3D"color:rgb(188= ,190,196)">/omega_c<br>x =3D irfft(X, </span><span style=3D"color:rgb(170,7= 3,38)">n</span><span style=3D"color:rgb(188,190,196)">=3Dn)<br>y =3D irfft(= </span><span style=3D"color:rgb(42,172,184)">2j</span><span style=3D"color:= rgb(188,190,196)">*np.pi*f*X, </span><span style=3D"color:rgb(170,73,38)">n= </span><span style=3D"color:rgb(188,190,196)">=3Dn)=C2=A0 </span><span styl= e=3D"color:rgb(122,126,133)"># dx / dt</span></p> <p>=C2=A0</p><p style=3D"margin:0px"># Two ShortTimeFFT instances are neede= d:</p> <p>=C2=A0<span style=3D"color:rgb(188,190,196)">kw =3D </span><span style= =3D"color:rgb(136,136,198)">dict</span><span style=3D"color:rgb(188,190,196= )">(</span><span style=3D"color:rgb(170,73,38)">hop</span><span style=3D"co= lor:rgb(188,190,196)">=3D</span><span style=3D"color:rgb(42,172,184)">10</s= pan><span style=3D"color:rgb(188,190,196)">, </span><span style=3D"color:rg= b(170,73,38)">fs</span><span style=3D"color:rgb(188,190,196)">=3D</span><sp= an style=3D"color:rgb(42,172,184)">1</span><span style=3D"color:rgb(188,190= ,196)">/T,</span><span style=3D"color:rgb(170,73,38)">fft_mode</span><span = style=3D"color:rgb(188,190,196)">=3D</span><span style=3D"color:rgb(106,171= ,115)">'onesided'</span><span style=3D"color:rgb(188,190,196)">, </= span><span style=3D"color:rgb(170,73,38)">phase_shift</span><span style=3D"= color:rgb(188,190,196)">=3D</span><span style=3D"color:rgb(207,142,109)">No= ne</span><span style=3D"color:rgb(188,190,196)">)<br>SFT =3D ShortTimeFFT.f= rom_dual(windows.hann(</span><span style=3D"color:rgb(42,172,184)">32</span= ><span style=3D"color:rgb(188,190,196)">, </span><span style=3D"color:rgb(1= 70,73,38)">sym</span><span style=3D"color:rgb(188,190,196)">=3D</span><span= style=3D"color:rgb(207,142,109)">False</span><span style=3D"color:rgb(188,= 190,196)">), **kw)</span></p> <p>=C2=A0<span style=3D"color:rgb(122,126,133)"># Differentiate dual window= per FFT:</span></p> <p>=C2=A0<span style=3D"color:rgb(188,190,196)">diff_dual_win =3D irfft(rff= t(SFT.dual_win) * </span><span style=3D"color:rgb(42,172,184)">2j</span><sp= an style=3D"color:rgb(188,190,196)">*np.pi*rfftfreq(SFT.m_num))<br>dSFT =3D= ShortTimeFFT.from_dual(diff_dual_win, **kw)</span></p> <p>=C2=A0</p><p>=C2=A0<span style=3D"color:rgb(122,126,133)"># Perform the = filtering:</span></p> <p>=C2=A0<span style=3D"color:rgb(188,190,196)">S_x =3D SFT.stft(x)<br>dS_x= =3D </span><span style=3D"color:rgb(42,172,184)">2j</span><span style=3D"c= olor:rgb(188,190,196)">*np.pi*SFT.f[:, np.newaxis] * S_x<br>dx =3D SFT.istf= t(dS_x, </span><span style=3D"color:rgb(170,73,38)">k1</span><span style=3D= "color:rgb(188,190,196)">=3Dn) + dSFT.istft(S_x, </span><span style=3D"colo= r:rgb(170,73,38)">k1</span><span style=3D"color:rgb(188,190,196)">=3Dn)</sp= an></p> <p>=C2=A0</p><p>=C2=A0<span style=3D"color:rgb(122,126,133)"># Plot windows= :</span></p> <p>=C2=A0<span style=3D"color:rgb(188,190,196)">fg0, axx0 =3D plt.subplots(= </span><span style=3D"color:rgb(42,172,184)">2</span><span style=3D"color:r= gb(188,190,196)">, </span><span style=3D"color:rgb(42,172,184)">1</span><sp= an style=3D"color:rgb(188,190,196)">, </span><span style=3D"color:rgb(170,7= 3,38)">sharex</span><span style=3D"color:rgb(188,190,196)">=3D</span><span = style=3D"color:rgb(106,171,115)">'all'</span><span style=3D"color:r= gb(188,190,196)">, </span><span style=3D"color:rgb(170,73,38)">tight_layout= </span><span style=3D"color:rgb(188,190,196)">=3D</span><span style=3D"colo= r:rgb(207,142,109)">True</span><span style=3D"color:rgb(188,190,196)">)<br>= axx0[</span><span style=3D"color:rgb(42,172,184)">0</span><span style=3D"co= lor:rgb(188,190,196)">].set(</span><span style=3D"color:rgb(170,73,38)">tit= le</span><span style=3D"color:rgb(188,190,196)">=3D</span><span style=3D"co= lor:rgb(106,171,115)">"Windows"</span><span style=3D"color:rgb(18= 8,190,196)">)<br>axx0[</span><span style=3D"color:rgb(42,172,184)">1</span>= <span style=3D"color:rgb(188,190,196)">].set(</span><span style=3D"color:rg= b(170,73,38)">title</span><span style=3D"color:rgb(188,190,196)">=3D</span>= <span style=3D"color:rgb(106,171,115)">"Dual Windows"</span><span= style=3D"color:rgb(188,190,196)">)<br>axx0[</span><span style=3D"color:rgb= (42,172,184)">0</span><span style=3D"color:rgb(188,190,196)">].plot(SFT.win= , </span><span style=3D"color:rgb(106,171,115)">'.-'</span><span st= yle=3D"color:rgb(188,190,196)">, </span><span style=3D"color:rgb(170,73,38)= ">alpha</span><span style=3D"color:rgb(188,190,196)">=3D</span><span style= =3D"color:rgb(42,172,184)">0.5</span><span style=3D"color:rgb(188,190,196)"= >, </span><span style=3D"color:rgb(170,73,38)">label</span><span style=3D"c= olor:rgb(188,190,196)">=3D</span><span style=3D"color:rgb(106,171,115)">= 9;SFT'</span><span style=3D"color:rgb(188,190,196)">)<br>axx0[</span><s= pan style=3D"color:rgb(42,172,184)">0</span><span style=3D"color:rgb(188,19= 0,196)">].plot(dSFT.win, </span><span style=3D"color:rgb(106,171,115)">'= ;.-'</span><span style=3D"color:rgb(188,190,196)">, </span><span style= =3D"color:rgb(170,73,38)">alpha</span><span style=3D"color:rgb(188,190,196)= ">=3D</span><span style=3D"color:rgb(42,172,184)">0.5</span><span style=3D"= color:rgb(188,190,196)">, </span><span style=3D"color:rgb(170,73,38)">label= </span><span style=3D"color:rgb(188,190,196)">=3D</span><span style=3D"colo= r:rgb(106,171,115)">'dSFT'</span><span style=3D"color:rgb(188,190,1= 96)">)<br>axx0[</span><span style=3D"color:rgb(42,172,184)">1</span><span s= tyle=3D"color:rgb(188,190,196)">].plot(SFT.dual_win, </span><span style=3D"= color:rgb(106,171,115)">'.-'</span><span style=3D"color:rgb(188,190= ,196)">, </span><span style=3D"color:rgb(170,73,38)">alpha</span><span styl= e=3D"color:rgb(188,190,196)">=3D</span><span style=3D"color:rgb(42,172,184)= ">0.5</span><span style=3D"color:rgb(188,190,196)">, </span><span style=3D"= color:rgb(170,73,38)">label</span><span style=3D"color:rgb(188,190,196)">= =3D</span><span style=3D"color:rgb(106,171,115)">'SFT'</span><span = style=3D"color:rgb(188,190,196)">)<br>axx0[</span><span style=3D"color:rgb(= 42,172,184)">1</span><span style=3D"color:rgb(188,190,196)">].plot(dSFT.dua= l_win, </span><span style=3D"color:rgb(106,171,115)">'.-'</span><sp= an style=3D"color:rgb(188,190,196)">, </span><span style=3D"color:rgb(170,7= 3,38)">alpha</span><span style=3D"color:rgb(188,190,196)">=3D</span><span s= tyle=3D"color:rgb(42,172,184)">0.5</span><span style=3D"color:rgb(188,190,1= 96)">, </span><span style=3D"color:rgb(170,73,38)">label</span><span style= =3D"color:rgb(188,190,196)">=3D</span><span style=3D"color:rgb(106,171,115)= ">'dSFT'</span><span style=3D"color:rgb(188,190,196)">)</span></p> <p>=C2=A0</p><p>=C2=A0<span style=3D"color:rgb(122,126,133)"># Plot signal:= </span></p> <p>=C2=A0<span style=3D"color:rgb(188,190,196)">fg1, ax1 =3D plt.subplots()= <br>ax1.set(</span><span style=3D"color:rgb(170,73,38)">title</span><span s= tyle=3D"color:rgb(188,190,196)">=3D</span><span style=3D"color:rgb(106,171,= 115)">rf"STFT-based Differentiator for $f_c=3D</span><span style=3D"co= lor:rgb(207,142,109)">{</span><span style=3D"color:rgb(188,190,196)">f[k_c]= </span><span style=3D"color:rgb(207,142,109)">}</span><span style=3D"color:= rgb(106,171,115)">\,$Hz Signal"</span><span style=3D"color:rgb(188,190= ,196)">,<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 </span><span style= =3D"color:rgb(170,73,38)">xlabel</span><span style=3D"color:rgb(188,190,196= )">=3D</span><span style=3D"color:rgb(106,171,115)">"Time $t$"</s= pan><span style=3D"color:rgb(188,190,196)">, </span><span style=3D"color:rg= b(170,73,38)">ylabel</span><span style=3D"color:rgb(188,190,196)">=3D</span= ><span style=3D"color:rgb(106,171,115)">"Amplitude"</span><span s= tyle=3D"color:rgb(188,190,196)">)<br>ax1.plot(t, omega_c * x, </span><span = style=3D"color:rgb(170,73,38)">alpha</span><span style=3D"color:rgb(188,190= ,196)">=3D</span><span style=3D"color:rgb(42,172,184)">0.5</span><span styl= e=3D"color:rgb(188,190,196)">, </span><span style=3D"color:rgb(170,73,38)">= label</span><span style=3D"color:rgb(188,190,196)">=3D</span><span style=3D= "color:rgb(106,171,115)">r"$x(t) \cdot 2\pi f_c$"</span><span sty= le=3D"color:rgb(188,190,196)">)<br>ax1.plot(t, y, </span><span style=3D"col= or:rgb(106,171,115)">'--'</span><span style=3D"color:rgb(188,190,19= 6)">, </span><span style=3D"color:rgb(170,73,38)">alpha</span><span style= =3D"color:rgb(188,190,196)">=3D</span><span style=3D"color:rgb(42,172,184)"= >0.5</span><span style=3D"color:rgb(188,190,196)">, </span><span style=3D"c= olor:rgb(170,73,38)">label</span><span style=3D"color:rgb(188,190,196)">=3D= </span><span style=3D"color:rgb(106,171,115)">"$y =3D dx/dt$"</sp= an><span style=3D"color:rgb(188,190,196)">)<br>ax1.plot(t, dx, </span><span= style=3D"color:rgb(170,73,38)">alpha</span><span style=3D"color:rgb(188,19= 0,196)">=3D</span><span style=3D"color:rgb(42,172,184)">0.5</span><span sty= le=3D"color:rgb(188,190,196)">, </span><span style=3D"color:rgb(170,73,38)"= >label</span><span style=3D"color:rgb(188,190,196)">=3D</span><span style= =3D"color:rgb(106,171,115)">"STFT-based $dx/dt$"</span><span styl= e=3D"color:rgb(188,190,196)">)</span></p> <p>=C2=A0</p><p>=C2=A0</p><p>=C2=A0<span style=3D"color:rgb(207,142,109)">f= or </span><span style=3D"color:rgb(188,190,196)">ax_ </span><span style=3D"= color:rgb(207,142,109)">in </span><span style=3D"color:rgb(188,190,196)">(*= axx0, ax1):<br>=C2=A0=C2=A0=C2=A0 ax_.legend()<br>=C2=A0=C2=A0=C2=A0 ax_.gr= id(</span><span style=3D"color:rgb(207,142,109)">True</span><span style=3D"= color:rgb(188,190,196)">)</span></p> <p>=C2=A0</p><p style=3D"margin:0px">plt.show()</p> <p>=C2=A0</p> <br><p style=3D"margin:0px">Cheers, dietrich</p> <br><br><br><p style=3D"margin:0px">[1] <a href=3D"https://ltfat.org/doc/ga= bor/" target=3D"_blank">https://ltfat.org/doc/gabor/</a></p> <p style=3D"margin:0px">[2] <a href=3D"https://scipy.github.io/devdocs/tuto= rial/signal.html#spectral-analysis" target=3D"_blank">https://scipy.github.= io/devdocs/tutorial/signal.html#spectral-analysis</a></p> <p style=3D"margin:0px">[3] <a href=3D"https://scipy.github.io/devdocs/tuto= rial/signal.html#short-time-fourier-transform" target=3D"_blank">https://sc= ipy.github.io/devdocs/tutorial/signal.html#short-time-fourier-transform</a>= </p> <br><p style=3D"margin:0px">On Saturday, 2 March 2024 22:10:45 CET Edward R= ichards wrote:</p> <p style=3D"margin:0px">> Thank you for your response, Dietrich.</p> <p style=3D"margin:0px">> </p> <p style=3D"margin:0px">> We are currently implementing the differentiat= ion as you suggested: (1)</p> <p style=3D"margin:0px">> STFT the signal, (2) multiply each slice by j2= =CF=80f, (3) ISTFT. We are not</p> <p style=3D"margin:0px">> using a rectangular window, though I can see h= ow it performs. I default to</p> <p style=3D"margin:0px">> a nuttall for DSP, and its relatively poor per= formance lead to my question.</p> <p style=3D"margin:0px">> Is there a way of predicting if an operation d= oes depend on the window? If</p> <p style=3D"margin:0px">> so, is there a standard practice for selection= ?</p> <p style=3D"margin:0px">> </p> <p style=3D"margin:0px">> This procedure does show promise. The integrat= ion step is much simpler with</p> <p style=3D"margin:0px">> the STFT than a simple FFT/IFFT.</p> <p style=3D"margin:0px">> </p> <p style=3D"margin:0px">> On Sat, Mar 2, 2024 at 1:49=E2=80=AFAM Dietric= h Brunn <<a href=3D"mailto:[email protected]" target=3D"_blank">Diet= [email protected]</a>> wrote:</p> <p style=3D"margin:0px">> > Hi Ned,</p> <p style=3D"margin:0px">> > </p> <p style=3D"margin:0px">> > it is nice to see that the new STFT funct= ionality is actually being used.</p> <p style=3D"margin:0px">> > </p> <p style=3D"margin:0px">> > > Our first use case is an integrator = and differentiator. We implement</p> <p style=3D"margin:0px">> > > this</p> <p style=3D"margin:0px">> > > </p> <p style=3D"margin:0px">> > > with a forward/inverse stft pair, an= d the frequency domain calculus</p> <p style=3D"margin:0px">> > > </p> <p style=3D"margin:0px">> > > definitions.</p> <p style=3D"margin:0px">> > </p> <p style=3D"margin:0px">> > May I ask how you implemented the differe= ntiation in the STFT space? Note</p> <p style=3D"margin:0px">> > that differentiation *is* window dependen= t. An ad-hoc approach would be</p> <p style=3D"margin:0px">> > using a rectangular dual window (utilizin= g `ShortTimeFFT.from_dual`) and</p> <p style=3D"margin:0px">> > multiplying the STFT by j2=CF=80f.</p> <p style=3D"margin:0px">> > </p> <p style=3D"margin:0px">> > Cheers, dietrich</p> <br><br><p></p><p></p><p></p><p></p><p></p><p></p><p></p><p></p></div> </blockquote></div> --000000000000d8cac40613027684-- --===============1865187239829479041== Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ SciPy-Dev mailing list -- [email protected] To unsubscribe send an email to [email protected] https://mail.python.org/mailman3/lists/scipy-dev.python.org/ Member address: [email protected] --===============1865187239829479041==--