bagustris@/home: Memfilter Sinyal Suara Pada Range Frekuensi Tertentu (Bandpass) dengan Matlab / Octave + Python

Sunday, March 30, 2014

Memfilter Sinyal Suara Pada Range Frekuensi Tertentu (Bandpass) dengan Matlab / Octave + Python

Buka - buka file lama di folder kuliah S1 Teknik Fisika, saya menemukan file script matlab yang saat itu dimintai tolong seorang teman untuk mengerjakan Tugas Akhirnya. Script ini untuk memfilter sinyal (suara) pada frekuensi tertentu dengan menggunakan bandpass filter, yakni gabungan antara low pass filter dan high-pass filter.

Misalnya kita mempunyai sinyal suara dengan frekuensi 0 - 20000 Hz, namun karena suatu hal (misalnya noise) kita cuma butuh sinyal di frekuensi 700 - 1400 Hz saja, maka kita bisa menggunakan script matlab / octave ini. Script ini berjalan lancar baik pada Matlab maupun GNU Octave.

Gambar spectrogram di bawah ini membandingkan spectrogram sinyal sebelum difilter dan setelah difilter. Script untuk GNU / Octave ada di bawahnya.

Sinyal Sebelum di filter

Sinyal Setelah di filter

Berikut script kode Octave / Matlab (link Github disini).

% demo program untuk memfilter suara pada range fekuensi oktaf tertentu
% by bagustris bagustris@linuxmail.org
% adopted from somewhere (forget aka lupa)

% set general variables
clear all; close all; clc;
sf = 22050;  % sample frequency
nf = sf / 2; % nyquist frequency
d = 1.0;     % duration, ganti untuk merubah durasi suara
n = sf * d;  % number of samples
nh = n / 2;  % half number of samples

% set variables for filter
lf = 700;   % lowest frequency, ganti dengan batas bawah bandpass
hf = 1400;   % highest frequency, ganti dengan batas atas bandpass
lp = lf * d; % ls point in frequency domain
hp = hf * d; % hf point in frequency domain

% band pass filter!
filter = zeros(1, n);           % initializaiton by 0
filter(1, lp : hp) = 1;         % filter design in real number
filter(1, n - hp : n - lp) = 1; % filter design in imaginary number

% make noise
rand('state',sum(100 * clock));  % initialize random seed
noise = randn(1, n);             % Gausian noise
noise = noise / max(abs(noise)); % -1 to 1 normalization

% do filter
s = fft(noise);                  % FFT
s = s .* filter;                 % filtering
s = ifft(s);                     % inverse FFT
s = real(s);

%untuk mendengarkan suara asli
sound(noise,sf);
pause(d+0.5);

%untuk mendengarkan suara terfilter (bandpass);
sound(s,sf);
pause(d+0.5);

%untuk  menyimpan file suara dalam format .wav
%ketik perintah berikut dalam command window (tanpa tanda persen%)
% wavwrite(s,sf 'namafile.wav');

Untuk mengecek hasilnya (spectrum yang telah terfilter) bisa menggunakan spectrogram sbb:

figure(1); specgram(noise, 512, sf);
figure(2); specgram(s, 512, sf);

Silakan dicoba. Please note the correctness and accurateness are not guaranteed.

Update 2020/06:
Dengan cara yang serupa, teknik yang sama dapat diterapkan dengan Python Scipy, yakni dengan filter butterworth. Berikut caranya (sumber: stackoverflow).

from scipy.signal import butter, lfilter

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y


def run():
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import freqz

    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 5000.0
    lowcut = 500.0
    highcut = 1250.0

    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [3, 6, 9]:
        b, a = butter_bandpass(lowcut, highcut, fs, order=order)
        w, h = freqz(b, a, worN=2000)
        plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)

    plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)],
             '--', label='sqrt(0.5)')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    # Filter a noisy signal.
    T = 0.05
    nsamples = int(T * fs)
    t = np.linspace(0, T, nsamples, endpoint=False)
    a = 0.02
    f0 = 600.0
    x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
    x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
    x += a * np.cos(2 * np.pi * f0 * t + .11)
    x += 0.03 * np.cos(2 * np.pi * 2000 * t)
    plt.figure(2)
    plt.clf()
    plt.plot(t, x, label='Noisy signal')

    y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)
    plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
    plt.xlabel('time (seconds)')
    plt.hlines([-a, a], 0, T, linestyles='--')
    plt.grid(True)
    plt.axis('tight')
    plt.legend(loc='upper left')

    plt.show()


run()

Hasilnya adalah sebagai berikut.