STFT và ISTFT nghịch đảo trong Python

Question

Có bất kỳ hình thức có mục đích chung nào là short-time Fourier transform với biến đổi nghịch đảo tương ứng được tích hợp vào SciPy hoặc NumPy hay bất kỳ thứ gì không?

Có những pyplot specgram chức năng trong matplotlib, trong đó kêu gọi ax.specgram(), trong đó kêu gọi mlab.specgram(), trong đó kêu gọi _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. 

nhưng

Đây là một chức năng helper mà thực hiện các sự tương đồng giữa 204 #psd, csd và phổ. Đó là KHÔNG nghĩa là để được sử dụng bên ngoài mLab

Tôi không chắc chắn nếu điều này có thể được sử dụng để làm một STFT và ISTFT, mặc dù. Có gì khác không, hoặc tôi nên dịch một cái gì đó như these MATLAB functions?

Tôi biết cách viết triển khai ad-hoc của riêng mình; Tôi chỉ đang tìm một cái gì đó đầy đủ tính năng, có thể xử lý các chức năng cửa sổ khác nhau (nhưng có một mặc định sane), hoàn toàn không thể đảo ngược với các cửa sổ COLA (istft(stft(x))==x), được thử nghiệm bởi nhiều người, không có lỗi từng cái một, xử lý kết thúc và không có đệm tốt, RFFT thực hiện nhanh chóng cho đầu vào sản, vv

Answer 1

Tôi là một chút muộn để điều này, nhưng nhận ra scipy có inbuilt istft chức năng như của 0.19.0

Answer 2

tôi cũng thấy điều này trên GitHub, nhưng có vẻ như để hoạt động trên đường ống thay vì mảng bình thường:

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

Đã tìm thấy một STFT khác, nhưng không có hàm nghịch đảo tương ứng:

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

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

• w là một chức năng cửa sổ như một mảng
• L là sự chồng chéo, trong các mẫu

Answer 4

Đây là mã Python của tôi, đơn giản hóa cho việc này trả lời:

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 

Ghi chú:

1. Danh sách hiểu là một mẹo nhỏ tôi thích sử dụng để mô phỏng xử lý khối tín hiệu trong NumPy/scipy. Nó giống như blkproc ở Matlab. Thay vì một vòng lặp for, tôi áp dụng một lệnh (ví dụ: fft) cho mỗi khung của tín hiệu bên trong một danh sách hiểu và sau đó scipy.array chuyển nó thành một mảng 2D. Tôi sử dụng điều này để tạo ra quang phổ, sắc ký, MFCC-gram, và nhiều hơn nữa.
2. Ví dụ này, tôi sử dụng phương pháp chồng chéo và thêm ngây thơ trong istft. Để tái tạo lại tín hiệu ban đầu, tổng của các hàm cửa sổ tuần tự phải là hằng số, tốt nhất là bằng sự thống nhất (1.0).Trong trường hợp này, tôi đã chọn cửa sổ Hann (hoặc hanning) và chồng chéo 50% hoạt động hoàn hảo. Xem this discussion để biết thêm thông tin.
3. Có thể có nhiều cách để tính toán ISTFT hơn. Ví dụ này chủ yếu có nghĩa là giáo dục.

Xét nghiệm:

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

Nếu bạn có quyền truy cập vào một thư viện nhị phân C mà những gì bạn muốn, sau đó sử dụng http://code.google.com/p/ctypesgen/ để tạo ra một giao diện Python đó thư viện.

Answer 6

Đây là mã STFT mà tôi sử dụng. STFT + ISTFT ở đây cung cấp cho tái thiết hoàn hảo (ngay cả đối với các khung hình đầu tiên). Tôi hơi sửa đổi mã được đưa ra ở đây bởi Steve Tjoa: ở đây độ lớn của tín hiệu tái tạo giống như tín hiệu đầu vào.

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

Không có câu trả lời nào ở trên hoạt động tốt OOTB cho tôi. Vì vậy, tôi đã sửa đổi Steve 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 và istft trông khá giống với những gì tôi đang tìm kiếm, mặc dù họ chưa hề tồn tại vào thời điểm đó:

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

Họ không đảo ngược chính xác, Tuy nhiên; các đầu được giảm dần.

Answer 9

Tôi nghĩ scipy.signal có những gì bạn đang tìm kiếm. Nó có giá trị mặc định hợp lý, hỗ trợ nhiều loại cửa sổ, vv ...

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)