Perfect reconstruction: Short-Time Fourier Transform (STFT) and Invers-nya menggunakan half-cycle sine window [2]

Tulisan kali ini sama persis dengan tulisan sebelumnya di link ini. Saya hanya menulis ulang dari Octave ke Python. Selamat mencoba STFT sederhana dengan Python.

 

Langkah pertama yakni membuat modul untuk short time fourier transform serta inversnya dengan nama "stft.py"

import numpy as np

def stft_sin(x,N,nfft):
 
  # Zero padding
  ncol = np.ceil(2*(x.size-N/2.)/N)
  if (ncol+1)*N/2. > x.size:
    temp = x
    x = np.zeros((ncol+1)*N/2)
    x[:temp.size] = temp
 
  temp = x
  x = np.zeros(temp.size+N)
  x[N/2:N/2+temp.size] = temp
 
  # Window sinus
  n = np.arange(N)
  Nindex = (x.size-N/2)*2/N
  window = np.tile(np.sin((n.reshape(n.size,1)+0.5)*np.pi/N),Nindex)
 
  # Proses windowing
  colindex = np.arange(Nindex)*N/2.
  rowindex = n.reshape(n.size,1)
  y = np.zeros((N,Nindex))
  y = x[(rowindex + colindex).astype(int)]
  y = y * window
 
  # FFT
  X = np.fft.fft(y,n=nfft,axis=0)
  X = X[:nfft/2+1]
 
  return X

def istft_sin(X,N,ndata):
  temp = X
  X = np.zeros((temp.shape[0]*2-2,temp.shape[1]),dtype=complex)
  X[:temp.shape[0],:temp.shape[1]] = temp
  X[temp.shape[0]::,:temp.shape[1]] = temp[np.arange(-2,-temp.shape[0],-1),:].conj()
  s = np.real(np.fft.ifft(X,axis=0))[:N,::]
  x = np.zeros((1,(s.shape[1]+1)*N/2))
  p = np.copy(x)
 
  # Window sinus
  n = np.arange(N)
  window = np.tile(np.sin((n.reshape(n.size,1)+0.5)*np.pi/N),(1,X.shape[1]))
 
  for i in range(X.shape[1]):
    x[0,i*N/2:(i+2)*N/2] = x[0,i*N/2:(i+2)*N/2] + (s[:,i]*window[:,i]).reshape(1,s.shape[0])
    p[0,i*N/2:(i+2)*N/2] = p[0,i*N/2:(i+2)*N/2] + (window[:,i]**2).reshape(1,s.shape[0])
 
  x = x/p
  x = x[0,N/2:ndata+N/2]
 
  return x

Untuk mengetahui apakah program kita sudah benar atau belum, kita akan membangkitkan sinyal berikut ini

$x(n) \begin{cases} \sin(2 \pi 1000 n) & \quad \text{if } \ \  t \geq 0 \ \ \& \ \ t < 1 \\ \sin(2 \pi 2000 n) & \quad \text{if } \ \ t \geq 1 \ \ \& \ \ t < 2 \\  \sin(2 \pi 3000 n) & \quad \text{if } \ \  t \geq 2 \ \ \& \ \ t < 3 \\   \end{cases}$  

dan mentransformasikannya ke domain waktu-frekuensi. Simpan file berikut ini dengan nama "demo_stft.py" dan jalankan di komputer anda. Hasilnya adalah seperti pada gambar di atas

import numpy as np
import matplotlib.pyplot as plt
import stft
import time

t0 = time.time()

fs = 10**4
t = np.arange(0,3,1./fs)
x = np.zeros(t.shape)

for i in range(3):
  range_t = (t >= i* t.max() /3.) & (t < (i+1)* t.max() /3.)
  x = x + np.sin(2*np.pi*(i+1)*1000*t) * range_t
 
# STFT
nfft = 1024
N = nfft/4.
X = stft.stft_sin(x,N,nfft)

#plot
plt.figure()
xy_label = [0,t.max(),0,fs/2.]
ratio = t.max()*2/fs
plt.imshow(np.abs(X),origin='lower',interpolation='nearest',cmap='jet',extent=xy_label,aspect=ratio)
plt.show()

t1 = time.time()
print 'Elapsed time %.5f seconds.' % (t1-t0)

Sama halnya dengan invers STFT, simpan file berikut dengan nama "demo_istft.py" dan jalankan lagi. Hasilnya seperti di bawah ini

import numpy as np
import matplotlib.pyplot as plt
import stft
import time

t0 = time.time()

fs = 10**4
t = np.arange(0,3,1./fs)
x = np.zeros(t.shape)
for i in range(3):
  range_t = (t >= i* t.max() /3.) & (t < (i+1)* t.max() /3.)
  x = x + np.sin(2*np.pi*(i+1)*5*t) * range_t
 
# STFT
nfft = 1024
N = nfft/4.
X = stft.stft_sin(x,N,nfft)

# ISTFT
ndata = t.size
x_new = stft.istft_sin(X,N,ndata)

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,x)
plt.ylim(-1.1,1.1)
plt.subplot(2,1,2)
plt.plot(t,x_new)
plt.ylim(-1.1,1.1)
plt.show()

print 'Error = %.5e' % np.max(np.abs(x-x_new))

t1 = time.time()
print 'Elapsed time %.5f' % (t1-t0)

Code untuk STFT dan inversnya dapat didownload di link ini.
Selamat mencoba