2. Sayısal Sistemlere Giriş

2.1. Sayısal Sistemlere Giriş

Önemli not: Bu defterin içeriği mühendislik öğrencileri/mezunlarına yönelik hazırlanmıştır. Okurların bahsi geçen konularda sinyal işleme kuramına aşina oldukları varsayılmıştır. Bu alanda deneyiminiz yok ise içerik sizin için yorucu olacaktır. Bu durumda 2. ve 3. defterleri atlayıp “Sinyallerin frekans spektrum analizi” başlıklı 4. defterden okumaya devam edebilirsiniz.

Bu defterimizde amacımız sayısal sinyalleri girdi olarak alıp işleyen ve sayısal çıktı sunan sistemlere örnekler sunmaktır. En basitinden başlayarak bazı sistemleri ele alacak, Python ile kodlayacak ve bu sistemlerin karakteristiklerini inceleyeceğiz.

Sistemler verilen bir girdiye karşılık bir çıktı üreten yapılar olarak ele alanabilir ve sıklıkla Şekilde 2.1’deki gibi bir kutu ile temsil edilir. Bu çizimde üstte sürekli zaman (analog işaret ve sistem), altta ayrık zaman (sayısal işaret ve sistem) temsil edilmiştir.

Şekil 2.1: Sistemler: üst grafik) sürekli zamanda bir sistem, alt grafik) ayrık zamanda bir (sayısal) sistem

Şekilde kutu içerisinde belirtilen \(h(t)\) (veya \(h[n]\)) sistemin bir dürtü(impulse) girdiye karşılık ürettiği çıktıyı temsil etmektedir.

\[h[n] = f(\delta[n])\]

Buna dürtü cevabı(impulse response) ismi verildiğini önceki defterimizde belirtmiştik.

Sistemlerimiz doğrusal ve zamanla değişmiyor (İng: linear time invariant, LTI) ve sistemin dürtü cevabı (\(h(t)\) veya \(h[n]\)) biliniyor ise herhangi bir girdi sinyaline karşılık sistemin oluşturacağı çıktı sinyali (\(y(t)\) veya \(y[n]\)) konvolüsyon işlemi ile hesaplanabilir (\(y[n] = h[n] \circledast x[n]\)). Bu işlemi de bir önceki defterimizde ele almıştık.

Sistemleri kod yazarak oluşturmak (bir bilgisayar programı olarak gerçeklemek) istediğimizde genellikle önümüzde dört olası yöntem bulunmaktadır:

• Tanımdan yola çıkarak hedeflenen sonucu temel vektör işlemleriyle veya döngüler yardımıyla gerçekleştiren bir fonksiyon yazmak

• Dürtü cevabını kullanarak konvolüsyon işlemi ile çıktıyı hesaplatmak

• Sistemin frekans cevabını polinomlarla (polinom katsayı dizileri) temsil etmek ve hazır sinyal işleme kütüphanesi fonksiyonlarını kullanmak (örnek: scipy.signal.lfilter)

• Sistemi frekans uzayında işlemlerle gerçeklemek

Burada öncelikle ilk iki yöntemi kullanarak bazı basit sistemler için Python fonksiyonları yazacağız. Daha sonra üçüncü ve dördüncü opsiyonlar için örnekler sunacağız.

Öncelikle gerekli kütüphanelerimizi yükleyerek başlayalım:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
%matplotlib inline
import IPython
from IPython.display import Image
from IPython.core.display import HTML 
import warnings
warnings.filterwarnings('ignore')

# Örneklerde kullanacağımız ses dosyasını indirelim
import urllib.request # Dosya indirmek için kullanacağımız kütüphane 
import soundfile as sf # Ses dosyalarını okumak için kullanacağımız kütüphane 

url = 'https://github.com/MTG/sms-tools/raw/master/sounds/speech-male.wav'
urllib.request.urlretrieve(url,'speech-male.wav');

2.1.1. Basit gecikme sistemi (Simple delay system)

Düşünebileceğimiz en basit sistemlerden birisi giriş sinyalini zamanda öteleyen/geciktiren bir sistemdir. Girdi-çıktı ilişkisi şu denklemle ifade edilebilir:

\[y[n] = f(x[n]) = x[n-n_0]\]

Bu sistem, \(x[n]\) girdi sinyalini \(n_0\) örnek kadar geciktirecektir. Sistemin dürtü cevabını bulmak istediğimizde \(\delta[n]\)in girdi olduğu durumda (diğer bir deyişle \(x[n]\) yerine \(\delta[n]\) yazdığımızda) oluşan çıktıyı hesaplayabiliriz:

\[h[n] = f(\delta[n]) = \delta[n-n_0]\]

Şimdi iki yöntem ile bu sayısal sistemi gerçekleyelim. Sistemi temsilen fonksiyon tanımlayacağız. \(x[n]\) girdi sinyali olarak fonksiyona verilecek ve fonksiyon sonuç olarak \(y[n]\) çıktı sinyalini döndürecek.

# Tanımdan yola çıkarak sistem gerçekleme
def geciktir_v1(x_n, n0, ayni_uzunlugu_kullan=False):
  '''Girdi sinyali olarak verilen x[n]'i n0 kadar öteleyerek
    çıktı oluşturup döndüren fonksiyon - sürüm 1'''
  # Sinyali kaydırmak için başına sıfırlar ekleyelim
  y_n = np.concatenate((np.zeros((n0,)), x_n), axis=None)
  
  # Çıktı boyunun belirlenmesi
  if ayni_uzunlugu_kullan:
    return y_n[:x_n.size]
  else:
    return y_n

# Dürtü cevabı kullanarak sistem gerçekleme
def geciktir_v2(x_n, n0, ayni_uzunlugu_kullan=False):
  '''Girdi sinyali olarak verilen x[n]'i n0 kadar öteleyerek
    çıktı oluşturup döndüren fonksiyon - sürüm 2'''
  # Dürtü cevabının oluşturulması
  N = n0 + 1 # dürtü cevabı uzunluğu
  h_n = np.zeros((N,))
  h_n[n0] = 1
  # Konvolüsyon ile hesaplama
  y_n = np.convolve(x_n, h_n)

  # Çıktı boyunun belirlenmesi
  if ayni_uzunlugu_kullan:
    return y_n[:x_n.size]
  else:
    return y_n

Örnek sinyal yükleyerek iki fonksiyonu test edelim

# Sinyalin dosyadan okunması
ses, ornekleme_fr = sf.read('speech-male.wav')

# Uzun ses sinyalinden elle seçilmiş bir kısmının alınması ve pencere fonksiyonu ile çarpılması
baslangic_endeksi = 8800; bitis_endeksi = 9500; # başlangıç ve bitiş noktası değerleri keyfi seçildi
x_n = ses[baslangic_endeksi:bitis_endeksi] * signal.get_window('hann', bitis_endeksi-baslangic_endeksi)

# Sinyali girdi olarak kullanarak geciktir_v1 ve geciktir_v2 fonksiyonlarının 
#  çıktılarını elde edip girdi sinyali ile beraber görselleştirelim
n0 = 100 # örnek sayısı cinsinden gecikme miktarı

fig = plt.figure(figsize=(8,6))

plt.subplot(2,1,1)
plt.plot(x_n, label='x[n]')
plt.plot(geciktir_v1(x_n, n0), label='y[n] (v1)') # geciktir_v1 çıktısı çizdiriliyor
plt.ylabel("geciktir_v1 çıktısı")
plt.legend();plt.grid();

plt.subplot(2,1,2)
plt.plot(x_n, label='x[n]')
plt.plot(geciktir_v2(x_n, n0), label='y[n] (v2)') # geciktir_v2 çıktısı çizdiriliyor
plt.ylabel("geciktir_v2 çıktısı")
plt.xlabel("Zaman (örnek sayısı/endeksi)")
plt.legend();plt.grid();

Şekil 2.2: Geciktirme sistemi girdi(\(x[n]\)) ve çıktı(\(y[n]\)) sinyalleri

İki fonksiyon ile de aynı çıktıyı aldık. geciktir_v1 ve geciktir_v2 fonksiyonları aynı sonucu üretiyor.

Sistemin dürtü cevabını biliyoruz:

\[h[n] = \delta[n-n_0]\]

Sistemin bir de frekans cevabını çizdirelim. Bunun için dürtü cevabının (\(h[n]\)) Fourier dönüşümünü numpy.fft.fft fonksiyonu ile elde edip çizdirelim.

# Numpy.fft.fft fonksiyonu kullanarak frekans cevabı çizdiren fonksiyon
def frekans_cevabi_cizdir(x_n, frekans_nokta_sayisi=512):
  # Fourier dönüşümü bize simetrik bir sonuç vereceği için bütününü oluşturduktan sonra ilk kısmını alalım
  # bu amaçla iki kat boyut kullanıp ortadan böleceğiz
  X_f = np.fft.fft(x_n, frekans_nokta_sayisi * 2)
  X_f = X_f[:frekans_nokta_sayisi]
  
  # Genlik ve faz bileşenlerinin hesaplanması
  genlik_w = np.abs(X_f)
  faz_w = np.angle(X_f)

  # Spektrum çizdirme adımları
  plt.figure(figsize=(12,3))

  plt.subplot(1,2,1)
  plt.title('Genlik spektrumu')
  plt.plot(genlik_w, 'b')
  plt.ylabel('|X(f)|', color='b')
  plt.xlabel('frekans [örnek endeksi]')
  plt.grid()
  
  plt.subplot(1,2,2)
  plt.title('Faz spektrumu (katlanmış)')
  plt.plot(faz_w, 'b')
  plt.ylabel('faz(X(f)) (radyan)', color='b')
  plt.xlabel('frekans [örnek endeksi]')
  plt.grid()
  plt.show()

\(h[n] = \delta[n-n_0]\) sinyalini 200 örnek içerecek şekilde oluşturup karşılık gelen frekans cevabını çizdirelim

n0 = 100 # örnek sayısı cinsinden gecikme miktarı
N = n0 * 2 # dürtü cevabı uzunluğu, deltayı içerecek şekilde yeterince uzun seçelim
h_n = np.zeros((N,))
h_n[n0] = 1
frekans_cevabi_cizdir(h_n)

Şekil 2.3: Geciktirme sistemi frekans cevabı

Çizdirilen frekans cevabını incelediğimizde gecikme sisteminin sinyalin her frekanstaki bileşenini 1 ile çarptığını (genliğini değiştirmediğini), sadece fazı değiştirdiğini görebiliyoruz. Faz oldukça karışık bir görünüme sahip. Bunun temel nedeni fazın \([-\pi, \pi]\) aralığında temsil edilmesi (diğer bir deyişle faz spektrumunun bu aralığa “katlanması” (warping)). Spektrumdaki gerçek bir faz değerinin \(6.5\pi\) olduğunu düşünün. Eldeki kompleks sayıdan fazı elde ederseniz sonucu \(0.5\pi\) olarak gözleyeceksiniz. İşte bu tüm spektruma uygulanmış, bu sebeple faz spektrumu \([-\pi, \pi]\) aralığına katlanmış olarak gösterilmektedir. Katlanmış spektrumu açmak çoğu zaman zor bir işlemdir ancak basit sinyaller için numpy.unwrap fonksiyonunu kullanmayı deneyebiliriz. Aşağıda frekans spektrum çizdirme fonksiyonumuzu bu adımı ekleyerek tekrar yazalım ve kullanalım:

# Numpy.fft.fft fonksiyonu kullanarak frekans cevabı çizdiren fonksiyon (faz spektrumunu numpy.unwrap ile açar)
def frekans_cevabi_cizdir(x_n, frekans_nokta_sayisi=512):
  # Fourier dönüşümü bize simetrik bir sonuç vereceği için bütününü oluşturduktan sonra ilk kısmını alalım
  # bu amaçla iki kat boyut kullanıp ortadan böleceğiz
  X_f = np.fft.fft(x_n, frekans_nokta_sayisi * 2)
  X_f = X_f[:frekans_nokta_sayisi]

  # Genlik ve faz bileşenlerinin hesaplanması
  genlik_w = np.abs(X_f)
  faz_w = np.unwrap(np.angle(X_f)) # "unwrap"(katlamayı açma) işlemi eklendi
  
  # Spektrum çizdirme adımları
  plt.figure(figsize=(12,3))
  
  plt.subplot(1,2,1)
  plt.title('Genlik spektrumu')
  plt.plot(genlik_w, 'b')
  plt.ylabel('|X(f)|', color='b')
  plt.xlabel('frekans [örnek endeksi]')
  plt.grid()
  
  plt.subplot(1,2,2)
  plt.title('Faz spektrumu')
  plt.plot(faz_w, 'b')
  plt.ylabel('faz(X(f)) (radyan)', color='b')
  plt.xlabel('frekans [örnek endeksi]')
  plt.grid()
  plt.show()

frekans_cevabi_cizdir(h_n)

Şekil 2.4: Geciktirme sistemi frekans cevabı (faz spektrumu katlaması açılmış (phase unwrapped))

Faz spektrumumuz doğrusallaştı. Faz spektrumunun doğrusal olması oldukça önemli bir özellik. Doğrusal fazlı bir sistem, girdi sinyalinin içerisindeki farklı frekansta bileşenlerin hepsinin zamanda aynı miktarda kaydıracaktır. Bu özelliği defterimizin sonuna doğru örnek üzerinden açıklıyoruz. Şimdi bir sonraki sistemimizi ele alalım.

2.1.2. Hareketli ortalama filtresi (Moving average filter)

Şimdi bir sinyalin küçük varyasyonlarını gözardı edip genel yönelimini izlememize yardımcı olan bir sistemi ele alalım. Finansal veriler üzerinde sıkça uygulanan bu filtre/sistem, verinin belirli aralıklarda ortalamalarını hesaplamaya dayanır. Örneğin elimizde saatlik Dolar/TL parite değerlerini içeren bir sinyal olduğunu düşünelim. Saatlik veriler çok fazla iniş-çıkış içereceği için günlük ortalama, haftalık ortalama, aylık ortalama alarak değişimleri gözlemek isteyebiliriz. Bu sistem girdi olarak aldığı sinyali belirlediğimiz genişlikte pencerelere bölüp, ortalama hesaplayıp, elde ettiği değerleri içeren yeni bir sinyali çıktı olarak verecektir. Sistemi şu formülle tanımlayabiliriz:

\[y[n] = \frac{1}{M_1+M_2+1} \sum_{k = -M_1}^{M_2} x[n-k]\]

Burada, çıktının \(n\). noktadaki değerini (\(y[n]\)) hesaplarken girdinin \(n\). noktadan \(M_1\) örnek öncesi ile \(M_2\) örnek sonrası arası değerlerinin ortalaması alınmaktadır. Şimdi ilk olarak sistemi gerçeklemek için, bu denklemdeki toplama işlemini gerçekleştiren bir fonksiyon yazalım:

# Tanımdan yola çıkarak sistem gerçekleme
def hareketli_ortalama_v1(x_n, M_1, M_2):
  y_n = np.zeros_like(x_n)
  # Sinyalin her bir noktasının çevresinde bir kesit alınıp ortalamasının hesaplanması
  for ind in range(M_1, x_n.size - M_2):
    y_n[ind] = np.sum(x_n[ind-M_1:ind+M_2+1]) / (M_1 + M_2 + 1)
  return y_n 

Dürtü cevabını oluşturarak da sistemimizi gerçekleyebiliyorduk, bunu da yapalım.

\[h[n] = \frac{1}{M_1+M_2+1} \sum_{k = -M_1}^{M_2} \delta[n-k]\]

\(h[n]\), \(n=-M_1\) ve \(n=+M_2\) arasında \(\delta\) sinyallerinin toplamında oluşuyor. Öncelikle \(h[n]\)i oluşturalım, çizdirelim ve inceleyelim:

N = 100 # sinyaldeki toplam örnek sayısı 
M_1 = 10; M_2 = 10; # rasgele seçilmiş M1 ve M2 değerleri
n = np.arange(- N // 2, N // 2)
h_n = np.zeros((N,))
h_n[ np.logical_and(-M_1 <= n , n <= M_2)] = 1
h_n = h_n / (M_1 + M_2 + 1)
fig = plt.figure(figsize=(8,3))
plt.title('Hareketli ortalama filtresi dürtü cevabı ($M_1$ = 10, $M_2$ = 10)')
plt.ylabel('h[n]');plt.xlabel('n [örnek endeksi]');
plt.stem(n, h_n);plt.grid();

Şekil 2.5: Hareketli ortalama sistemi dürtü cevabı

Şimdi sistemimizi dürtü cevabı oluşturup girdi sinyaliyle konvolüsyon gerçekleştirecek şekilde yazalım:

# Dürtü cevabı kullanarak sistem gerçekleme
def hareketli_ortalama_v2(x_n, M_1, M_2):
  # Dürtü cevabının oluşturulması
  N = M_1 + M_2 + 1 # sinyaldeki toplam örnek sayısı 
  h_n = np.ones((N,))
  h_n = h_n / N

  # Çıktıyı konvolüsyon ile hesaplama
  y_n = np.convolve(x_n, h_n)
  
  # Konvolüsyon sonucu sinyal uzadı, keserek x_n ile aynı uzunlukta döndürelim
  return y_n[M_1:M_1+x_n.size]

Şimdi yazdığımız iki fonksiyonla elde ettiğimiz çıktıları karşılaştıralım:

# Sinyali girdi olarak kullanarak hareketli_ortalama_v1 ve hareketli_ortalama_v2 
#. fonksiyonlarının çıktılarını elde edip girdi sinyali ile beraber görselleştirelim
M_1 = 10; M_2 = 10;

fig = plt.figure(figsize=(8,6))

plt.subplot(2,1,1)
plt.plot(x_n, label='x[n]')
plt.plot(hareketli_ortalama_v1(x_n, M_1, M_2), label='y[n] (v1)') # hareketli_ortalama_v1 çıktısı çizdiriliyor
plt.ylabel("hareketli_ortalama_v1 çıktısı")
plt.legend();plt.grid();

plt.subplot(2,1,2)
plt.plot(x_n, label='x[n]')
plt.plot(hareketli_ortalama_v2(x_n, M_1, M_2), label='y[n] (v2)') # hareketli_ortalama_v2 çıktısı çizdiriliyor
plt.ylabel("hareketli_ortalama_v2 çıktısı")
plt.xlabel("Zaman (örnek sayısı/endeksi)")
plt.legend();plt.grid();

Şekil 2.6: Hareketli ortalama sistemi girdi(\(x[n]\)) ve çıktı(\(y[n]\)) sinyalleri

İki fonksiyonla aynı sonucu elde ettiğimizi doğrulamış olduk. Şimdi sistemin frekans cevabını inceleyelim:

frekans_cevabi_cizdir(h_n)

Şekil 2.7: Hareketli ortalama sistemi frekans cevabı

Sistemimizin düşük geçiren bir filtre olduğunu görebiliyoruz (sistem düşük frekans bileşenlerini yüksek frekans bileşenlerine görece daha yüksek sayılarla çarpıyor). Bir sonraki sistem örneğimize geçelim.

2.1.3. Kümülatif toplam yapan sistem (Accumulator function)

Basit sistemlerden bir diğeri çıktı sinyalinin \(n\) anındaki değerini hesaplamak için girdi sinyalinin \(n\) ve \(n\)’den küçük (geçmiş) değerlerinin tümünü toplayan sistemdir. Bu, bir çeşit integral olarak da düşünülebilir.

\[y[n] = \sum_{k = -\infty}^{n} x[k]\]

Sistemin dürtü cevabını bulmak için her zaman olduğu gibi girdi sinyali yerine delta fonksiyonunu yazabiliriz:

\[x[n] = \delta[n]\]

Sonuçta elde edilen dürtü cevabı:

\[h[n] = u[n]\]

olarak hesaplanabilir.

Sistemimizin dürtü cevabı birim basamakmış. Birim basamak sinyalini ilk defterimizin başında oluşturmuştuk. Aynı örneği tekrar inceleyelim:

\[\begin{split}u[n] = \begin{cases} 0, & \quad n < 0\\ 1, & \quad n \geq 0 \end{cases} \end{split}\]

N = 100 # seçilen örnek sayısı
n = np.arange(- N // 2, N // 2)
h_n_integral = np.zeros((N, )) 
h_n_integral[n >= 0] = 1

# Oluşturduğumuz dürtü cevabını çizdirelim
fig = plt.figure(figsize=(12,3))
plt.title('Sistem dürtü cevabı')
plt.ylabel('h[n]');plt.xlabel('n [örnek endeksi]');
plt.stem(n, h_n_integral);plt.grid();

Şekil 2.8: Akümülatör sistemi dürtü cevabı

Frekans cevabını da çizdirelim:

frekans_cevabi_cizdir(h_n_integral)

Şekil 2.9: Akümülatör sistemi frekans cevabı

Yine düşük geçiren bir sistem elde ettik. Burada önemli bir nokta sistemin sonsuz uzunlukta bir dürtü cevabı (IIR: Infinite impulse response) olması ancak bizim \(h[n]\) vektörünü sınırlı uzunlukta tutmak zorunda olmamızdır. Zamanda sınırlandırma yaptığımız durumda sistemin karakterini ve spektrumunu da değiştiririz. IIR sistemleri bu sebeple dürtü cevabı oluşturarak gerçeklememeyi tercih ederiz. Alternatif gerçekleme yöntemini ileride ele alacağız. Bu sistem özelinde: bu sistem stabil bir sistem olmadığı (dürtü cevabı toplamı sonlu olmadığı) için zaten sistemin frekans cevabını (Fourier dönüşümünü) hesaplayıp çizdiremeyiz (bakınız)). Yalnızca dürtü cevabını sınırlandırdığımız durumda dönüşüm hesaplayıp bir çizim elde edebiliyoruz.

2.1.4. Fark/türev sistemi (Differentiator)

Yukarıda integral alan sistem tanımlamıştık. Şimdi ters işlem olan fark alan sistemi tanımlayalım:

\[y[n] = x[n] - x[n-1]\]

Bu sistem \(n\) anındaki çıktıyı, girdi sinyalinin o andaki değerinin bir önceki değerinden farkı olarak hesaplar. Dürtü cevabını hesaplamak için \(x[n] = \delta[n]\) yazarsak:

\[h[n] = \delta[n] - \delta[n-1]\]

Dürtü cevabını oluşturalım ve spektrumuyla beraber çizdirelim:

N = 30
n = np.arange(- N // 2, N // 2)
h_n_fark = np.zeros((N, )); 
h_n_fark[n==0] = 1; h_n_fark[n==1] = -1 
fig = plt.figure(figsize=(12,3))
plt.stem(n, h_n_fark);plt.grid();plt.title('Fark sistemi dürtü cevabı');
plt.ylabel('h[n]');plt.xlabel('n [örnek endeksi]');

Şekil 2.10: Fark sistemi dürtü cevabı

Bu sistemin dürtü cevabı sınırlı uzunluğa sahip (FIR: finite impulse response). Sıfırları dışarıda bırakarak temsil etmemizde bir sakınca yok.

h_n_fark = np.array([1,-1])
frekans_cevabi_cizdir(h_n_fark)

Şekil 2.11: Fark sistemi frekans cevabı

Fark sisteminin yüksek geçiren bir filtre karakteristiği olduğunu gözlemlemiş olduk. İntegral alan sistem ile tam ters karakteristiğe sahip olmasını beklerken gözlemimiz bunu doğrulamadı. Bunun sebebi IIR olan integral sistemini sınırlı uzunluktaki bir dürtü cevabı ile temsil etmemiz oldu.

Verilen iki sistemin birbirinin tersi olup olmadığını seri bağlandıklarında birbirlerinin etkisini yok edip etmediklerini test ederek inceleyebiliriz. Ardışık bağlanan bu iki sistemi tek bir sistem olarak ele alıp, bu birleşik sistemin aldığı girdi sinyalini bozmadan çıktıda verip vermediğini kontrol edilebiliriz. (şekil değiştirilecek)

Şekil 2.12: Ters sistemlerin seri bağlanması ile elde edilen sistem

Şunu kontrol etmek istiyoruz:

\(y[n] = h^{-1}[n]\circledast(h[n] \circledast x[n]) = (h^{-1}[n]\circledast h[n]) \circledast x[n] =? x[n]\)

Bu doğruluğun sağlanabilmesi için;

\(h^{-1}[n] \circledast h[n] = \delta[n]\) olmalı. Bunu elimizdeki akümülatör ve fark sistemleri için kontrol edelim;

\(h^{-1}[n] \circledast h[n] = u[n] \circledast (\delta[n] - \delta[n-1])\) \( = u[n] \circledast \delta[n] - u[n] \circledast \delta[n-1]\)

\( = u[n] - u[n-1]\)

\( = \delta[n]\)

Evet, bu iki sistemin birbirinin tersi olduğunu doğruladık. Sinyal örneği üzerinden de bir test yapabiliriz. İntegral sistemini sadece pozitif zamandaki değerlerinden ibaret oluşturursak elimizdeki sinyal üzerinde bunu aşağıdaki gibi doğrulamayı deneyebiliriz.

\(x[n]\) → integal → fark → ? \(x[n]\)

Girdi ve çıktı sinyallerini beraber çizdirelim ve ne düzeyde örtüştüklerini görsel olarak incelemeye çalışalım. Çizimde üstüste bineceği için sinyallerden birisini y ekseninde 0.1 birim kaydırarak çizdireceğiz. (siz bu bileşeni dışarıda bırakıp deneyi tekrarlayabilirsiniz)

h_n_integral = np.ones_like(x_n)
h_n_fark = np.array([1,-1])
# x[n] -> integral 
ilk_sistem_y_n = np.convolve(x_n, h_n_integral)
# -> fark -> ...
ikinci_sistem_y_n = np.convolve(ilk_sistem_y_n, h_n_fark)

fig = plt.figure(figsize=(8,3))
plt.plot(x_n+0.1, label='girdi sinyali + 0.1') # üstüste çakışmayı önleme için sinyale 0.1 eklendi
plt.plot(ikinci_sistem_y_n[:x_n.size], label='ardışık olarak integral ve fark uygulanmış sinyal')
plt.legend();plt.grid();

Şekil 2.13: Birbirinin tersi olan sistemlerin seri bağlanması sonucu elde edilen sistemin girdi ve çıktı sinyalleri. (girdi sinyali y ekseninde 0.1 kaydırılmıştır)

Evet, sistemler seri bağlanarak oluşturulan sistemin girdi ile aynı çıktıyı ürettiğini bu örnekte görmüş olduk.

İki sistemin dürtü cevaplarının konvolüsyonunu hesaplayıp delta fonksiyon elde edip etmediğimizi de kontrol edebilirdik:

N = 50
n = np.arange(- N // 2, N // 2)

h_n_integral = np.zeros((N, )); h_n_integral[n >= 0] = 1
h_n_fark = np.zeros((N, )); h_n_fark[0] = 1;h_n_fark[1] = -1

sonuc = np.convolve(h_n_integral, h_n_fark)
sonuc = sonuc[:N]

fig = plt.figure(figsize=(12,3))
plt.stem(n, sonuc);plt.title('İntegral ve fark sistemlerinin seri bağlanması sonucu oluşan sistemin dürtü cevabı');

Şekil 2.14: Bir sistemin dürtü cevabının ters sistemin dürtü cevabı ile konvolüsyonu sonucu elde edilen bileşik sistem dürtü cevabı

Bu testte de, iki sistemin birbirinin tersi olduğu beklentimiz (sonuç \(\delta[n]\) çıktığı için) doğrulandı.

2.2. Egzersizler/Örnekler

2.2.1. Egzersiz/örnek 1

Girdi sinyalini 10 ile çarparak yükselten bir sistemin frekans cevabını çizdiren bir kod parçası yazınız.

\[y[n] = 10x[n]\]

\[h[n] = 10\delta[n]\]

carpan = 10
imp_resp_duration_n = 10
delta_n = np.zeros((imp_resp_duration_n,))
delta_n[0] = 1
h_n = carpan * delta_n
frekans_cevabi_cizdir(h_n)

Şekil 2.15: Sabit kazanç uygulayan sistem frekans cevabı

Sistem tüm frekanslardaki bileşenlerin genliklerini 10 ile çarpıyor ve fazlarını değiştirmiyor.

2.2.2. Egzersiz/örnek 2

Nedensel bir sistemin dürtü cevabı aşağıdaki gibi verilmiştir:

\(h[n]\) = [0.02455383, 0.23438946, 0.4821134, 0.23438946, 0.02455383]

Sistemin frekans cevabını çizdiriniz. Sistem düşük geçiren mi yüksek geçiren midir?

h_n = np.array([0.02455383, 0.23438946, 0.4821134, 0.23438946, 0.02455383])
fig = plt.figure(figsize=(6,3))
plt.stem(h_n);plt.grid();plt.title('Dürtü cevabı')
frekans_cevabi_cizdir(h_n)

Şekil 2.16: Egzersiz 2’de dürtü cevabı verilen sistemin frekans cevabı

Sistem düşük geçiren bir filtre karakteristiği göstermektedir. Bu örneği hazırlamak için dürtü cevabını hazır bir filtre tasarım fonksiyonu kullanarak oluşturduk. Altta kullandığımız kodu bulabilirsiniz:

uzunluk = 5 # h[n]'in uzunluğu
f_kesim = 0.25 # kesim frekansı frekans bandının ilk çeyreğinde seçilmiş
h_n = signal.firwin(uzunluk, f_kesim)
print(h_n)

[0.02455383 0.23438946 0.4821134  0.23438946 0.02455383]

Yukarıdaki örnekteki filtrenin hazır bir FIR filtre tasarım fonksiyonu kullanılarak oluşurulduğunu görmüş olduk. Siz de ihtiyaç duyduğunuzda signal.firwin fonksiyonu ile arzu ettiğiniz kesim frekansına sahip filtre tasarlayıp filtrenin dürtü cevabını elde edebilirsiniz. Burada önemli detaylardan birisi de oluştural filtrenin faz spektrumunun doğrusal olması; uygulandığı sinyalin değişik frekanslardaki bileşenlerini zamanda eşit miktarda kaydırmasıdır. Tasarlanan sistemin dürtü cevabının zamanda simetrik olması bu özelliği garanti altına alır. Bu konuyu defterin ilerleyen bölümlerinde ele alıyoruz.

2.2.3. Egzersiz/örnek 3

Örnek olarak sırasıyla bir düşük geçiren, bir yüksek geçiren, bir de bant geçiren filtre tasarlayıp elimizdeki örnek ses sinyaline uygulayıp sonucu dinleyip sistemlerin etkilerini deneyimleyelim.

# Konuşma sinyaline 3000 Hz kesim frekanslı düşük geçiren filtre uygulanması
uzunluk = 101 # h[n]'in uzunluğu
f_max = ornekleme_fr //2 # Örnekleme kuramına göre sinyalde temsil edilen maksimum frekans
f_kesim = 3000 # kesim frekansını 3000 Hz olarak belirleyelim
f_norm = f_kesim / f_max # scipy.signal.firwin "cutoff" parametresi, maksimum frekansa oran olarak 0-1 aralığında belirleniyor
LP_h_n = signal.firwin(uzunluk, f_norm, pass_zero = "lowpass")
frekans_cevabi_cizdir(LP_h_n)
filtrelenmis_ses = np.convolve(ses, LP_h_n)
IPython.display.Audio(filtrelenmis_ses, rate=ornekleme_fr)

Şekil 2.17: Düşük geçiren filtre frekans cevabı

Yukarıdaki üçgen “play” butonuna tıklayarak filtre çıktısı elde edilen sesi dinleyebilirsiniz.

# Konuşma sinyaline 1000 Hz kesim frekanslı yüksek geçiren filtre uygulanması
uzunluk = 101
f_max = ornekleme_fr //2 # Örnekleme kuramına göre sinyalde temsil edilen maksimum frekans
f_kesim = 1000 # kesim frekansını 1000 Hz olarak belirleyelim
f_norm = f_kesim / f_max
LP_h_n = signal.firwin(uzunluk, f_norm, pass_zero = "highpass")
frekans_cevabi_cizdir(LP_h_n)
filtrelenmis_ses = np.convolve(ses, LP_h_n)
IPython.display.Audio(filtrelenmis_ses, rate=ornekleme_fr)