FIRフィルタをPythonとＣ言語で実装する

概要

信号処理にノイズ除去という視点でフィルタがよく使われる。その中でも一番シンプルなFIRフィルタをPythonとC言語で実装する。Pythonベースでフィルタの評価を行なった上で、その結果をC言語ベースのアプリケーションに組み込む流れを想定する。

Pythonコード

以下にscipyを使用したFIRフィルタの実装を示す。元信号xは乱数とし、scipy.signal.firwinによってフィルタ係数を求める。scipyによるフィルタ適用後の信号y及び、べた打ちによるフィルタ適用後の信号dを比較する。

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

import numpy as np
import scipy.signal
from pylab import *
import csv
import matplotlib.pyplot as plt

#xはなんでもいい
x = np.random.rand(512)
plt.plot(x)
plt.show()

fs = 512 #サンプリング周波数
nyq = fs / 2.0  #ナイキスト周波数

# フィルタの設計
f1 = 1 / nyq      #カットオフ周波数1
f2 = 30.0 / nyq      #カットオフ周波数2
numtaps = 255       #フィルタ係数の数(奇数)

b = scipy.signal.firwin(numtaps, f1)                         # ローパス
#b = scipy.signal.firwin(numtaps, f2, pass_zero=False)        # ハイパス
#b = scipy.signal.firwin(numtaps, [f1, f2], pass_zero=False) # バンドパス

#フィルタ係数の保存
with open('coeff.csv', 'w') as f:
    writer = csv.writer(f, lineterminator='\n') # 改行コード（\n）を指定しておく
    writer.writerow(b)

# scipyによるFIRフィルタ
y = scipy.signal.lfilter(b, 1, x)     

# ベタ書きによるFIRフィルタ
d = np.zeros(fs)

for i in range(len(x)):
    d[i] = 0
    for j in range(len(b)):
        if(i-j)>=0:
            d[i] += b[j]*x[i-j]

plt.plot(y,color='r')
plt.plot(d,color='b')
plt.show()

元信号

フィルタ適用後（赤:y, 青:d）

Cコード

Python上でのべた打ち結果が正しければ、C言語コードに落とし込む。少々雑だが、bには、Pythonコードで生成したcoeff.csvの配列を格納する。csv形式で保存してあるので、テキストエディタ等で開けば、そのままコピー&ペーストで貼り付けられる。

    int sampleRate = 512

    x= // 雑だがｘを適当に記述する

    float* filteredData = (float*)malloc(sizeof(float)*512);

    float b[] = {} //雑だが、Pythonで作成したcoeff.csvの中身を括弧内に挿入する

    for(int i=0;i<sampleRate;i++){
        float d = 0;
        for(int j=0; j<(sizeof b)/(sizeof b[0]); j++){
            if((i-j)>=0){
                d += b[j]*x[i-j];
            }
        }
        filteredData[i]=d;
    }

    free(filteredData)

非常に単純だが備忘録程度なのでお許しを、FIRフィルタってどうやって実装するんだっけ？という人の助けになれば。