デジタルフィルタ（ローパスフィルタ）をめちゃくちゃ簡単に実装しよう（C言語） – FIRフィルタ編

こんにちは。Yukiです。めちゃくちゃ久しぶりのブログ投稿です。

思ったより簡単にデジタルフィルタを実装できたので、備忘録として残しておこうと思います。
今回は、PC（C言語）を使い、音声ファイルをを使います。

注意：ここに書いてあることが正しいのかは知りません！参考程度にどうぞ！

めちゃ簡実装 – 1回目
デジタルフィルタの種類
今回使うツール
そもそもフィルタの設計とは何をするのか
実際に設計してみる
C言語で使えるファイル形式に変換
矩形波を生成する
デジタルフィルタを実装してみよう
波形から見るFIRフィルタの欠点
実際に音声データでやってみる
感想

めちゃ簡実装 – 1回目

何をやるかは未定ですが、今のところ次のような形で進める予定です。

1回目 – FIRフィルタを実装
- この記事でやります
2回目 – IIRフィルタを実装
- https://denshi1996.com/?p=843
3回目 – マイコンに実装
- 準備中…
4回目 – DSP機能付きマイコンで実装
- 準備中…

デジタルフィルタの種類

おそらくいろいろあるとは思うのですが、大きく分けて2つあります。

FIRフィルタ (Finite Impulse Response フィルタ)
- 常に安定
- 遅延が常に一定
- 設計がIIRフィルタに比べて簡単
IIRフィルタ (Infinite Impulse Response フィルタ)
- フィードバックがある
  - 不安定になる可能性がある
- 遅延が短い
- アナログフィルタと共通点が多い

簡単にいうとこんな感じです。今回はFIRフィルタを使います。（IIRフィルタは次回以降にやってみたいと思います）

あと、詳しいパラメータの計算方法とかはやりません（ツールを使います）
デジタルフィルタについて詳しく知りたい方は、別の記事をどうぞ

今回使うツール

デジタルフィルタアナライザ – DIGITALFILTER.COM
- https://digitalfilter.com/products/dfalz/jpdfalz.html

MATLAB等に入っているフィルタ設計ツールでも良いですが、無料のため今回はこれを使います。
株式会社デジタルフィルター様に感謝いたします。

VSCode

ぶっちゃけエディタならなんでいいです。今回はVSCodeを使用します。（デバッグ作業も楽だしね）

Audacity

波形を見るときに使用します。RAW形式の音声を取り込めれば何でも大丈夫です。

今回はGCCを利用します。（まあ多分C言語をコンパイルできるならなんでも大丈夫）

そもそもフィルタの設計とは何をするのか

自分はこれがさっぱりわからなかったです。結論からいうと、各重みの決める事がフィルタの設計になります。

FIRフィルタはこのような構成図になっています。

Zとかいうブロックがありますが、これは遅延ブロックで、前のデータを表しています。下にいけばいくほど、古いデータが流れていくみたいなイメージです。

FIRフィルタでは、h0とかh1とかh2とかの数値（パラメータ）を調整するというのが目標になります。この値を変えることによって、指定されたローパス・ハイパス・バンドパスフィルタなどを作ることができます。

（株価とかの平均線なども実は構成図自体は同じです）

実際に設計してみる

まずは、フィルタ設計ツールにて、係数（重み）を計算します。

Filter TypeをLPFにして、FIRにセットして、Designをクリックします。

今回は、次のように設定しました。

今回は1000Hzベースの矩形波を入力とするので、1000Hzちょっと大きいくらいの値をfc1(Pass band edge)にします。fc2(Stop band edge)は2000Hz以下くらいにします。(矩形波には理論上ですが3ω～の周波数成分が入っているので)

Tap Countは、増やすと繰り返し回数が増える（遅延が増大する）ので程々にします。（当たり前ですが、大きくしたほうがフィルタの性能は良くなりますが）

グラフ部分にフィルタ性能が表示されます。今回だと、カット帯域は-50dB以下 (大体0.003倍くらい)まで軽減されます。

次に、この特性になる係数を出力されます。

Coefficientsをクリック。

Quantize byにチェックを入れます。

Save/Loadをクリック。

一番上のSave Coefficientsにチェックを入れて、OKを押して、どっかに保存してください。（名前はなんでも大丈夫です）

C言語で使えるファイル形式に変換

ここからVSCodeを利用します。

実はMATLAB等のツールを使うと、直接C言語用のヘッダファイルが生成されるのですが、このツールはそのような機能がないため、ヘッダファイルに改造します。

VSCodeに取り込み、名前を〇〇◯.h (C++で作る場合は◯◯◯.hpp)に変更します。

今回は FIR_Filter_LPF.h という名前にしました。

Ctrl + H を押すと置換機能が出ます。　上のテキストボックスには置換前の文字列。下のテキストボックスには置換後の文字列を入れます。

ここがVSCodeの便利機能なのですが、■*みたいなマークのところを押すと、\n等の改行文字も置換対象にすることができます。

今回は、すべての行末にカンマを付与したいので、置換前に\n 置換後に ,\n を入れます。

#define宣言でTAP_Nを宣言します。（#defineを使えない派閥の方は、const int 等で適当に代用してください。そっちのほうがいいらしいですし）

次に、配列を作ります。 const double filter_weight[]という配列名にしました。

（画像のようにしてくれれば大丈夫）

最終行のカンマを1つ消し、 }; をつけて終了です。（セミコロンを忘れないように）

矩形波を生成する

矩形波を出力するためのプログラムを作ります。

#include <stdio.h>
#include <stdint.h>

//保存先
#define FILE_PATH "C:/Users/User/Desktop/square_wave.raw"

//サンプリング周波数
#define Fs 44100
//矩形波の周波数
#define SQ_F 1000
//矩形波を何秒生成するか
#define SQ_TIME 1
//矩形波のMAX値（符号あり16bitだったら32767）
#define SQ_MAX 25000
//矩形波のMIN値（符号あり16bitだったら-32768）
#define SQ_MIN -25000

int main(void){
    FILE *file;
    file = fopen(FILE_PATH, "wb");

    //矩形波を生成するために必要なサンプリング数を計算
    int samples_per_cycle = (int)(Fs / SQ_F);

    for(int j = 0; j < SQ_F * SQ_TIME ; j++){
        //半周期分の矩形波を生成
        for(int i = 0; i < samples_per_cycle / 2;i++){
            int16_t data = SQ_MAX; //HIGH
            fwrite(&data, sizeof(int16_t), 1, file);
        }

        for(int i = 0; i < samples_per_cycle / 2;i++){
            int16_t data = SQ_MIN; //LOW
            fwrite(&data, sizeof(int16_t), 1, file);
        }
    }

    fclose(file);
}

5行目の#define宣言は、任意の保存先に変更してください。

次に、この波形データをAudacityに取り込みます。（ちゃんと矩形波が生成されているのかの確認）

ファイル → インポート → Raw データをインポートをクリックします。

上に貼り付けたソースコードをそのまま使い、RAWデータを生成した場合はこのような設定になります。

こんな感じで矩形波が取り込めていたらOKです。（1秒に達していませんが、これはサンプリング周波数が44.1kHzで中途半端なため、丸め誤差が発生しているからです）

デジタルフィルタを実装してみよう

フィルタの係数・入力データの用意ができました。ということでデジタルフィルタをちゃちゃっと実装しましょう。

#include <stdio.h>
#include <stdint.h>
#include "FIR_Filter_LPF.h"

#define INPUT_FILE_PATH "C:/Users/User/Desktop/square_wave.raw"
#define OUTPUT_FILE_PATH "C:/Users/User/Desktop/filter_output.raw"

#define VOLUME 0.7

int main(void){
    FILE *input_file;
    FILE *output_file;

    //フィルタ用に過去の入力データを保存しておく配列
    int16_t input_buffer[TAP_N]= {0};

    input_file = fopen(INPUT_FILE_PATH, "rb");

    if(input_file == NULL){
        printf("INPUTファイルが開けませんでした\n");
        return -1;
    }

    output_file = fopen(OUTPUT_FILE_PATH, "wb");

    if(output_file == NULL){
        printf("OUTPUTファイルが開けませんでした\n");
        return -2;
    }

    while(1){
        int16_t input_data;
        size_t read_size = fread(&input_data, sizeof(int16_t), 1, input_file);

        if(read_size != 1){
            //読み込み失敗 (EOF(ファイルの最後)に到達したと考えて終了)
            break;
        }

        double output_data = 0.0;

        //バッファをずらす
        for(int i = TAP_N - 1; i > 0; i--){
            input_buffer[i] = input_buffer[i - 1];
        }
        input_buffer[0] = input_data * VOLUME;

        //畳み込み演算
        for(int i = 0; i < TAP_N; i++){
            output_data += filter_weight[i] * input_buffer[i];
        }

        //出力データのクリップ処理
        if(output_data > 32767.0){
            output_data = 32767.0;
        } else if(output_data < -32768.0){
            output_data = -32768.0;
        }

        int16_t output_data_int16 = (int16_t)output_data;

        fwrite(&output_data_int16, sizeof(int16_t), 1, output_file);
    }


}

実はフィルタの実装自体は下に示す部分のみです。

        //バッファをずらす
        for(int i = TAP_N - 1; i > 0; i--){
            input_buffer[i] = input_buffer[i - 1];
        }
        input_buffer[0] = input_data * VOLUME;

        //畳み込み演算
        for(int i = 0; i < TAP_N; i++){
            output_data += filter_weight[i] * input_buffer[i];
        }

思ったよりもシンプルじゃないですか？

このプログラムを実行すると、フィルタ後のRAWデータが出力されます。先ほどと同様の方法でAudacityに読み込ませます。