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

こんにちは。Yukiです。今回はIIRフィルタを実装してみたいと思います。

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

めちゃ簡フィルタ実装 – 2回目
今回使うツール
IIRフィルタの構造
実際に設計してみる
Cで使えるヘッダファイルに変換
プログラム
Audacityで確認
IIRフィルタも思ったより簡単だった

めちゃ簡フィルタ実装 – 2回目

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

今回使うツール

前回と変わりませんが、ここに出しておきます。

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

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

VSCode

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

Audacity

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

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

IIRフィルタの構造

IIRフィルタには様々な構造があります。

今回使用するツールが出力できるものはDirect Form I型という最も理論的にわかりやすいIIRフィルタです。

次数を増やすには、フィルタ出力の横に同じ形のブロック図をつける感じです。

（上の写真で構成するとより性能が良いフィルタになるが、計算量も増える）

今回のソフトでは最大8までいけるようです。

FIRフィルタでは200とかにしましたが、IIRフィルタはフィードバックもあり効率の良いフィルタなので、8も入れればかなり急峻なフィルタになります。

（IIRフィルタはInfinite Impulse Response フィルタという通り、無限の応答なので、フィードバックで前の情報を残して無限にしているらしいです）

実際に設計してみる

では早速dfalzを開き、設計してみましょう。

LPF・IIRにセットして、Designボタンを押します。

IIRフィルタの層数と、リプル、カットオフ周波数とサンプリング周波数を設定します。

今回は次のような感じで設定しました。（1000Hzの矩形波を入力するので、少し大きめにしています）

Coefficientsを選択します。

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

Save/Loadを押し、Save coefficientsをクリック。OKを押して任意の場所に保存すればOKです。

Cで使えるヘッダファイルに変換

FIRフィルタのときは、後ろにカンマを付けるだけだったので、VSCodeの置換機能を使えば簡単でしたが、IIRフィルタは様々な定数があるので若干面倒です。

ということで、このソフト形式で出てきたテキストファイル→ヘッダファイルになるCプログラムを作ります。

#include <stdio.h>
#include <string.h>

#define INPUT_FILE_PATH "IIR_LPF.txt"
#define OUTPUT_FILE_PATH "output.h"

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

    char buffer[256];

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

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

    fprintf(output_file, "#define IIR_A 0\n#define IIR_B 1\ndouble IIR_Weight[2][8][3];\nvoid IIR_Weight_setup(void)\n{\n");

    fgets(buffer, sizeof(buffer), input_file); //1行飛ばす

    int filter_cnt = 0;
    int AorB_cnt = 0;
    int filter_n = 0;

    while(fgets(buffer, sizeof(buffer), input_file) != NULL){

        //IIR_AまたはIIR_Bに直す
        char tmp[6];
        if(AorB_cnt == 0){
            snprintf(tmp,sizeof(tmp),"IIR_A");
        }else{
            snprintf(tmp,sizeof(tmp),"IIR_B");
        }

        //改行を消す
        buffer[strlen(buffer) - 1] = '\0';
        fprintf(output_file, "\tIIR_Weight[%d][%s][%d] = %s;\n", filter_cnt, tmp, filter_n, buffer);

        filter_n++;
        if(filter_n >= 3 && AorB_cnt == 0){
            filter_n = 1;
            AorB_cnt = 1;
        }
        if(filter_n >= 3 && AorB_cnt == 1){
            filter_n = 0;
            AorB_cnt = 0;
            filter_cnt++;
        }
    }

    fprintf(output_file, "}\n");

    fclose(input_file);
    fclose(output_file);
}

ちゃんと作っていないので、バグとかあるかもしれません。自己責任でお使いください。

4, 5行目にあるdefine宣言のところを変更するだけで使えます。

INPUT_FILE_PATHには、先ほどdfalz1で出力したテキストファイルのパスを入力します。OUTPUT_PATHは保存先を入力します。（上書きで保存されるので注意！）

プログラム

FIRフィルタのときのプログラムでも使いましたが、矩形波を出力するCプログラムはそのまま再利用します。

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

//保存先
#define FILE_PATH "filter_test.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);
}

上のdefine宣言のところが設定です。適宜設定して使ってください。（といっても保存先ぐらいしか触るところない気もしますが）

それでは、IIRフィルタのプログラムを作ります。といってもブロック図通りにCプログラムを書けばよいだけです。

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

#define INPUT_FILE_PATH "music.raw"
#define OUTPUT_FILE_PATH "music_output.raw"

#define IIR_N 8

#define VOLUME 0.7

#define BUFFER_A 0
#define BUFFER_B 1

double IIR_buffer[2][IIR_N][3] = {0};

int main(void)
{
    IIR_Weight_setup();

    FILE *input_file;
    FILE *output_file;

    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;
        }

        for (int i = 0; i < IIR_N; i++)
        {

            double output_data = 0.0;

            // バッファをずらす
            IIR_buffer[BUFFER_A][i][2] = IIR_buffer[BUFFER_A][i][1];
            IIR_buffer[BUFFER_A][i][1] = IIR_buffer[BUFFER_A][i][0];
            if (i == 0)
            {
                IIR_buffer[BUFFER_A][i][0] = input_data * VOLUME;
            }
            else
            {
                IIR_buffer[BUFFER_A][i][0] = IIR_buffer[BUFFER_B][i - 1][0];
            }

            IIR_buffer[BUFFER_B][i][2] = IIR_buffer[BUFFER_B][i][1];
            IIR_buffer[BUFFER_B][i][1] = IIR_buffer[BUFFER_B][i][0];

            output_data = IIR_Weight[IIR_A][i][0] * IIR_buffer[BUFFER_A][i][0];
            output_data += IIR_Weight[IIR_A][i][1] * IIR_buffer[BUFFER_A][i][1];
            output_data += IIR_Weight[IIR_A][i][2] * IIR_buffer[BUFFER_A][i][2];
            output_data += IIR_Weight[IIR_B][i][1] * IIR_buffer[BUFFER_B][i][1] * (-1);
            output_data += IIR_Weight[IIR_B][i][2] * IIR_buffer[BUFFER_B][i][2] * (-1);

            if (IIR_buffer[BUFFER_B][IIR_N - 1][0] > 32767.0)
            {
                IIR_buffer[BUFFER_B][IIR_N - 1][0] = 32767.0;
            }

            if (IIR_buffer[BUFFER_B][IIR_N - 1][0] < -32768.0)
            {
                IIR_buffer[BUFFER_B][IIR_N - 1][0] = -32768.0;
            }

            IIR_buffer[BUFFER_B][i][0] = output_data;
        }

        int16_t output_int16 = (int16_t)IIR_buffer[BUFFER_B][IIR_N - 1][0];
        fwrite(&output_int16, sizeof(int16_t), 1, output_file);
    }
}