Inwertowalne STFT i ISTFT w Pythonie

Question

Czy istnieje jakakolwiek uniwersalna forma short-time Fourier transform z odpowiednią transformacją odwrotną wbudowaną w SciPy lub NumPy czy coś?

Tam jest pyplot specgram funkcja w matplotlib, który nazywa ax.specgram(), który nazywa mlab.specgram(), który nazywa _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. 

ale

Jest to funkcja pomocnika, który implementuje wspólność między 204 #psd, csd i spektrogramem. Jest NIE przeznaczone do pracy na zewnątrz mLAB

Nie jestem pewien, czy to może być używane do zrobić STFT i ISTFT, choć. Czy jest coś jeszcze, czy powinienem przetłumaczyć coś w stylu these MATLAB functions?

Wiem, jak napisać własną implementację ad-hoc; Po prostu szukam czegoś w pełni wyposażonego, które może obsłużyć różne funkcje okienkowania (ale ma normalne domyślne), jest w pełni odwracalne z oknami COLA (istft(stft(x))==x), testowane przez wiele osób, bez błędów "jeden po drugim", obsługuje kończy i zero wyściółka dobrze, szybka realizacja RFFT za prawdziwe wejście itp

Answer 1

Jestem trochę późno na to, ale sobie sprawę scipy posiada wbudowane istft funkcję od 0.19.0

Answer 2

znalazłem też to na GitHub, ale wydaje się działać na rurociągach zamiast normalnych tablic:

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)) 

Answer 3

Znaleziono inne STFT, ale bez odpowiadającej funkcji odwrotnej:

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

def stft(x, w, L=None): 
    ... 
    return X_stft 

• w jest funkcją okienkującą tablicy
• L jest nakładanie w próbkach

Answer 4

Tutaj jest kod Python uproszczony tego odpowiedź:

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 

Uwagi:

1. listowego jest mały trick lubię używać do symulacji przetwarzania sygnałów w bloku numpy/scipy. To jest jak blkproc w Matlabie. Zamiast pętli for, stosuję polecenie (np. fft) do każdej ramki sygnału wewnątrz rozumienia listy, a następnie scipy.array rzuca ją do tablicy 2D. Używam tego do tworzenia spektrogramów, chromagramów, MFCC-gramów i wielu innych.
2. Dla tego przykładu używam naiwnej metody overlap-and-add w istft. Aby zrekonstruować oryginalny sygnał, suma sekwencyjnych funkcji okna musi być stała, najlepiej równa jedności (1,0).W tym przypadku wybrałem okno Hann (lub hanning) i 50% zachodzenie, które działa idealnie. Aby uzyskać więcej informacji, patrz this discussion.
3. Istnieją prawdopodobnie bardziej zasadnicze sposoby obliczania ISTFT. Ten przykład ma być głównie edukacyjny.

Test:

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)') 

Answer 5

Jeśli masz dostęp do C binarnym biblioteki nie to, co chcesz, a następnie użyć http://code.google.com/p/ctypesgen/ wygenerować interfejs Pythona do tego biblioteka.

Answer 6

Oto kod STFT, którego używam. STFT + ISTFT daje tutaj idealną rekonstrukcję (nawet dla pierwszych klatek). Lekko zmodyfikowałem kod podany tutaj przez Steve'a Tjoa: tutaj wielkość rekonstruowanego sygnału jest taka sama jak sygnału wejściowego.

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 

Answer 7

Żadna z powyższych odpowiedzi nie sprawdziła się dla mnie. Więc zmodyfikowałem Steve'a Tjoa.

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... 

Answer 8

librosa.core.stft i istft wyglądają bardzo podobnie do tego, co szukałem, choć nie istniały w chwili:

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

Nie odwracać dokładnie, chociaż; końce są zwężone.

Answer 9

Myślę, że scipy.signal ma to, czego szukasz. To ma uzasadnione domyślne, obsługuje wiele typów okien, itp ...

http://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.signal.spectrogram.html

from scipy.signal import spectrogram 
freq, time, Spec = spectrogram(signal)