Librosa 的 STFT 功能齐全,因此除非您非常小心地处理频谱,否则您不会从它的istft
.
这是我从头开始编写的一对函数stft
和,它们表示正向和反向 STFT,以及一个帮助方法,可为您提供 STFT 阵列中每个像素的时间和频率位置,以及一个演示:istft
import numpy as np
import numpy.fft as fft
def stft(x, Nwin, Nfft=None):
"""
Short-time Fourier transform: convert a 1D vector to a 2D array
The short-time Fourier transform (STFT) breaks a long vector into disjoint
chunks (no overlap) and runs an FFT (Fast Fourier Transform) on each chunk.
The resulting 2D array can
Parameters
----------
x : array_like
Input signal (expected to be real)
Nwin : int
Length of each window (chunk of the signal). Should be ≪ `len(x)`.
Nfft : int, optional
Zero-pad each chunk to this length before FFT. Should be ≥ `Nwin`,
(usually with small prime factors, for fastest FFT). Default: `Nwin`.
Returns
-------
out : complex ndarray
`len(x) // Nwin` by `Nfft` complex array representing the STFT of `x`.
See also
--------
istft : inverse function (convert a STFT array back to a data vector)
stftbins : time and frequency bins corresponding to `out`
"""
Nfft = Nfft or Nwin
Nwindows = x.size // Nwin
# reshape into array `Nwin` wide, and as tall as possible. This is
# optimized for C-order (row-major) layouts.
arr = np.reshape(x[:Nwindows * Nwin], (-1, Nwin))
stft = fft.rfft(arr, Nfft)
return stft
def stftbins(x, Nwin, Nfft=None, d=1.0):
"""
Time and frequency bins corresponding to short-time Fourier transform.
Call this with the same arguments as `stft`, plus one extra argument: `d`
sample spacing, to get the time and frequency axes that the output of
`stft` correspond to.
Parameters
----------
x : array_like
same as `stft`
Nwin : int
same as `stft`
Nfft : int, optional
same as `stft`
d : float, optional
Sample spacing of `x` (or 1 / sample frequency), units of seconds.
Default: 1.0.
Returns
-------
t : ndarray
Array of length `len(x) // Nwin`, in units of seconds, corresponding to
the first dimension (height) of the output of `stft`.
f : ndarray
Array of length `Nfft`, in units of Hertz, corresponding to the second
dimension (width) of the output of `stft`.
"""
Nfft = Nfft or Nwin
Nwindows = x.size // Nwin
t = np.arange(Nwindows) * (Nwin * d)
f = fft.rfftfreq(Nfft, d)
return t, f
def istft(stftArr, Nwin):
"""
Inverse short-time Fourier transform (ISTFT)
Given an array representing the output of `stft`, convert it back to the
original samples.
Parameters
----------
stftArr : ndarray
Output of `stft` (or something the same size)
Nwin : int
Same input as `stft`: length of each chunk that the STFT was calculated
over.
Returns
-------
y : ndarray
Data samples corresponding to STFT data.
See also:
stft : the forward transform
"""
arr = fft.irfft(stftArr)[:, :Nwin]
return np.reshape(arr, -1)
if __name__ == '__main__':
sampleRate = 100.0 # Hertz
N = 1024
Nwin = 64
# Generate a chirp: start frequency at 5 Hz and going down at 2 Hz/s
time = np.arange(N) / sampleRate # seconds
x = np.cos(2 * np.pi * time * (5 - 2 * 0.5 * time))
# Test with Nfft bigger than Nwin
Nfft = Nwin * 2
s = stft(x, Nwin, Nfft=Nfft)
y = istft(s, Nwin)
# Make sure the stft and istft are inverses. Caveat: `x` and `y` won't be
# the same length if `N/Nwin` isn't integral!
maxerr = np.max(np.abs(x - y))
assert (maxerr < np.spacing(1) * 10)
# Test `stftbins`
t, f = stftbins(x, Nwin, Nfft=Nfft, d=1 / sampleRate)
assert (len(t) == s.shape[0])
assert (len(f) == s.shape[1])
try:
import pylab as plt
plt.imshow(np.abs(s), aspect="auto", extent=[f[0], f[-1], t[-1], t[0]])
plt.xlabel('frequency (Hertz)')
plt.ylabel('time (seconds (start of chunk))')
plt.title('STFT with chirp example')
plt.show()
except ModuleNotFoundError:
pass
如果这对您来说更容易阅读,那么这是一个要点。
整个模块假定只有真实数据并使用 Numpy 的rfft
函数。您绝对可以将其推广到复杂数据(或使用 librosa),但对于您的应用程序(音频屏蔽),使用纯实数转换可以更轻松地确保一切正常,并且逆 STFT 的输出是纯实数(如果您正在执行完全通用的复杂 STFT,这很容易搞砸,您需要小心保持对称性)。
演示首先生成一些测试数据,并确认数据的istft
onstft
再次生成数据。测试数据是一个以 5 Hz 开始并以每秒 2 Hz 的速度下降的啁啾,因此在大约 10 秒的数据中,啁啾的频率会环绕并以 15 Hz 左右结束。该演示绘制了 STFT(通过获取 STFT 数组的绝对值):

所以
- 将此代码放入
stft.py
文件中,
- 将其导入为
import stft
,
- 将 STFT 计算为
spectrum = stft.stft(y, 128)
,
- 如演示中所示可视化您的频谱(确保添加
stft.
到定义的函数stft.py
!),
- 选择您想要衰减/放大的频率并将这些效果应用到
spectrum
阵列上,之前
- 最终通过
back_y = stft.istft(spectrum, 128)
.
掩蔽/放大/衰减频率内容意味着仅缩放spectrum
阵列的一些箱。如果您对如何执行此操作有具体问题,请告诉我们。但这有望为您提供一种应用任意效果的万无一失的方法。
如果您真的想使用 librosa 的功能,请告诉我们,我们也可以向您展示如何做到这一点。