#概要
信号処理にノイズ除去という視点でフィルタがよく使われる。その中でも一番シンプルな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= // 雑だが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フィルタってどうやって実装するんだっけ?という人の助けになれば。