Обратимый STFT и ISTFT в Python

Существует ли какая-либо универсальная форма кратковременного преобразования Фурье с соответствующим обратным преобразованием, встроенным в SciPy или NumPy или что-то еще?

Здесь функция pyplot specgram в matplotlib, которая вызывает ax.specgram(), которая вызывает mlab.specgram(), которая вызывает _spectral_helper():

#The checks for if y is x are so that we can use the same function to
#implement the core of psd(), csd(), and spectrogram() without doing
#extra calculations.  We return the unaveraged Pxy, freqs, and t.

но

Это вспомогательная функция, которая реализует общность между 204 #psd, csd и спектрограмму. это НЕ предназначен для использования вне mlab

Я не уверен, что это можно использовать для STFT и ISTFT. Есть ли что-нибудь еще, или я должен перевести что-то вроде эти функции MATLAB?

Я знаю, как написать свою собственную специальную реализацию; Я просто ищу что-то полнофункциональное, которое может обрабатывать различные функции окон (но имеет нормальное значение по умолчанию), полностью обратимо с окнами COLA (istft(stft(x))==x), проверенными несколькими людьми, без ошибок, хорошо обрабатывает концы и нулевое заполнение, быструю реализацию RFFT для реального ввода и т.д.

Ответы

Ответ 1

Я немного опаздываю на это, но реализованный scipy имеет встроенную istft с 0.19.0

Ответ 2

Вот мой код Python, упрощенный для этого ответа:

import scipy, pylab

def stft(x, fs, framesz, hop):
    framesamp = int(framesz*fs)
    hopsamp = int(hop*fs)
    w = scipy.hanning(framesamp)
    X = scipy.array([scipy.fft(w*x[i:i+framesamp]) 
                     for i in range(0, len(x)-framesamp, hopsamp)])
    return X

def istft(X, fs, T, hop):
    x = scipy.zeros(T*fs)
    framesamp = X.shape[1]
    hopsamp = int(hop*fs)
    for n,i in enumerate(range(0, len(x)-framesamp, hopsamp)):
        x[i:i+framesamp] += scipy.real(scipy.ifft(X[n]))
    return x

Примечания:

  • Понимание списка - это небольшой трюк, который мне нравится использовать для имитации блочной обработки сигналов в numpy/scipy. Это похоже на blkproc в Matlab. Вместо цикла for я применяю команду (например, fft) к каждому кадру сигнала внутри понимания списка, а затем scipy.array передает его в 2D-массив. Я использую это, чтобы создавать спектрограммы, хромаграммы, MFCC-граммы и многое другое.
  • В этом примере я использую метод наивного перекрытия и добавления в istft. Чтобы восстановить исходный сигнал, сумма последовательных оконных функций должна быть постоянной, предпочтительно равной единице (1.0). В этом случае я выбрал окно Hann (или hanning) и перекрытие 50%, которое отлично работает. Подробнее см. обсуждение.
  • Есть, вероятно, более принципиальные способы вычисления ISTFT. Этот пример в основном предназначен для обучения.

Тест:

if __name__ == '__main__':
    f0 = 440         # Compute the STFT of a 440 Hz sinusoid
    fs = 8000        # sampled at 8 kHz
    T = 5            # lasting 5 seconds
    framesz = 0.050  # with a frame size of 50 milliseconds
    hop = 0.025      # and hop size of 25 milliseconds.

    # Create test signal and STFT.
    t = scipy.linspace(0, T, T*fs, endpoint=False)
    x = scipy.sin(2*scipy.pi*f0*t)
    X = stft(x, fs, framesz, hop)

    # Plot the magnitude spectrogram.
    pylab.figure()
    pylab.imshow(scipy.absolute(X.T), origin='lower', aspect='auto',
                 interpolation='nearest')
    pylab.xlabel('Time')
    pylab.ylabel('Frequency')
    pylab.show()

    # Compute the ISTFT.
    xhat = istft(X, fs, T, hop)

    # Plot the input and output signals over 0.1 seconds.
    T1 = int(0.1*fs)

    pylab.figure()
    pylab.plot(t[:T1], x[:T1], t[:T1], xhat[:T1])
    pylab.xlabel('Time (seconds)')

    pylab.figure()
    pylab.plot(t[-T1:], x[-T1:], t[-T1:], xhat[-T1:])
    pylab.xlabel('Time (seconds)')

STFT of 440 Hz sinusoidISTFT of beginning of 440 Hz sinusoidISTFT of end of 440 Hz sinusoid

Ответ 3

Вот код STFT, который я использую. STFT + ISTFT дает идеальную реконструкцию (даже для первых кадров). Я слегка изменил код, приведенный здесь Стивом Тхоа: здесь величина восстановленного сигнала такая же, как и входного сигнала.

import scipy, numpy as np

def stft(x, fftsize=1024, overlap=4):   
    hop = fftsize / overlap
    w = scipy.hanning(fftsize+1)[:-1]      # better reconstruction with this trick +1)[:-1]  
    return np.array([np.fft.rfft(w*x[i:i+fftsize]) for i in range(0, len(x)-fftsize, hop)])

def istft(X, overlap=4):   
    fftsize=(X.shape[1]-1)*2
    hop = fftsize / overlap
    w = scipy.hanning(fftsize+1)[:-1]
    x = scipy.zeros(X.shape[0]*hop)
    wsum = scipy.zeros(X.shape[0]*hop) 
    for n,i in enumerate(range(0, len(x)-fftsize, hop)): 
        x[i:i+fftsize] += scipy.real(np.fft.irfft(X[n])) * w   # overlap-add
        wsum[i:i+fftsize] += w ** 2.
    pos = wsum != 0
    x[pos] /= wsum[pos]
    return x

Ответ 4

librosa.core.stft и istft выглядят довольно похожими на то, что я искал, хотя их не было в то время:

librosa.core.stft(y, n_fft=2048, hop_length=None, win_length=None, window=None, center=True, dtype=<type 'numpy.complex64'>)

Они не инвертируют точно; концы сужаются.

Ответ 5

Нашел еще один STFT, но не имеет соответствующей обратной функции:

http://code.google.com/p/pytfd/source/browse/trunk/pytfd/stft.py

def stft(x, w, L=None):
    ...
    return X_stft
  • w - это оконная функция в виде массива
  • L - это перекрытие, в образцах

Ответ 6

Ни один из вышеперечисленных ответов не работал хорошо для OOTB для меня. Поэтому я изменил Стив Тьоа.

import scipy, pylab
import numpy as np

def stft(x, fs, framesz, hop):
    """
     x - signal
     fs - sample rate
     framesz - frame size
     hop - hop size (frame size = overlap + hop size)
    """
    framesamp = int(framesz*fs)
    hopsamp = int(hop*fs)
    w = scipy.hamming(framesamp)
    X = scipy.array([scipy.fft(w*x[i:i+framesamp]) 
                     for i in range(0, len(x)-framesamp, hopsamp)])
    return X

def istft(X, fs, T, hop):
    """ T - signal length """
    length = T*fs
    x = scipy.zeros(T*fs)
    framesamp = X.shape[1]
    hopsamp = int(hop*fs)
    for n,i in enumerate(range(0, len(x)-framesamp, hopsamp)):
        x[i:i+framesamp] += scipy.real(scipy.ifft(X[n]))
    # calculate the inverse envelope to scale results at the ends.
    env = scipy.zeros(T*fs)
    w = scipy.hamming(framesamp)
    for i in range(0, len(x)-framesamp, hopsamp):
        env[i:i+framesamp] += w
    env[-(length%hopsamp):] += w[-(length%hopsamp):]
    env = np.maximum(env, .01)
    return x/env # right side is still a little messed up...

Ответ 7

Я также нашел это на GitHub, но, похоже, работает на конвейерах вместо обычных массивов:

http://github.com/ronw/frontend/blob/master/basic.py#LID281

def STFT(nfft, nwin=None, nhop=None, winfun=np.hanning):
    ...
    return dataprocessor.Pipeline(Framer(nwin, nhop), Window(winfun),
                                  RFFT(nfft))


def ISTFT(nfft, nwin=None, nhop=None, winfun=np.hanning):
    ...
    return dataprocessor.Pipeline(IRFFT(nfft), Window(winfun),
                                  OverlapAdd(nwin, nhop))

Ответ 9

Если у вас есть доступ к бинарной библиотеке C, которая делает то, что вы хотите, используйте http://code.google.com/p/ctypesgen/ для создания интерфейса Python для этой библиотеки.