firwin_test.py

#!/usr/bin/env python2
# -*- coding: utf-8 -*-

import numpy as np

from scipy.signal import firwin, firwin2, freqz
from scipy.io import loadmat

import matplotlib.pyplot as plt

b_matlab = loadmat('/home/larsoner/matlab/ir.mat')['b'][0]

fs = 1000.
n = len(b_matlab)
# Design with spectral inversion
# (e.g., http://mpastell.com/2010/01/18/fir-with-scipy/)
# This works better than "pass_zero=False" high-pass mode.
b_firwin = -firwin(n, cutoff=0.5 / (fs), window='hamming')
b_firwin[n // 2] = b_firwin[n // 2] + 1

nyq = fs / 2.
b_firwin2 = firwin2(n, [0, 0.5, nyq], [0, 1, 1], window='hamming', nyq=nyq)

fig, axes = plt.subplots(2)
for filt, label in zip((b_matlab, b_firwin, b_firwin2),
                       ('EEGLAB', 'firwin', 'firwin2')):
    w, h = freqz(filt, worN=10000)
    w *= fs / (2 * np.pi)
    axes[0].plot(w, np.abs(h), label=label)
    axes[1].plot(w, 20 * np.log10(np.abs(h)), label=label)
axes[0].legend()
axes[0].axhline(0.5, color='k', linestyle=':')
axes[1].axhline(-6, color='k', linestyle=':')
for ax in axes:
    ax.set(xlim=(0, 1))
axes[0].set(ylim=(0, 1), ylabel='Amplitude')
axes[1].set(ylim=(-100, 0), ylabel='Amplitude (dB)')
fig.tight_layout()