An IIR filter.

iir.erl

-module(iir).
-export([gen_filter/2, filter/3]).
-export([loop/4, aggregate_filtered/2, apply_filter/3]).
-compile({nowarn_unused_function, [aggregate/1,
                                   aggregate_buffered_acc/3,
                                   aggregate_buffered/2,
                                   aggregate_filtered_acc/3]}).
-define(PI, math:pi()).
-include_lib("eunit/include/eunit.hrl").

loop(DelayerA, DelayerB, PastX, PastY) ->
    receive
        {From, stop} ->
            From ! {self(), stop};
        {From, X} ->
            PastX2 = queue:liat(queue:cons(X, PastX)),
            Y = lists:sum(lists:map(fun({I, J}) -> I * J end,
                                    lists:zip(DelayerB, queue:to_list((PastX2))))),
            Y2 = Y + lists:sum(lists:map(fun({I, J}) -> -1 * I * J end,
                                         lists:zip(DelayerA, queue:to_list(PastY)))),
            From ! {self(), Y2},
            PastY2 = queue:liat(queue:cons(Y2, PastY)),
            loop(DelayerA, DelayerB, PastX2, PastY2)
    end.

start_loop(DelayerA, DelayerB) ->
    PastX = queue:from_list(lists:duplicate(length(DelayerB), 0)),
    PastY = queue:from_list(lists:duplicate(length(DelayerA), 0)),
    Pid = spawn_link(?MODULE, loop, [DelayerA, DelayerB, PastX, PastY]),
    Pid.

aggregate(Pid) ->
    receive
        {Pid, stop} ->
            ok;
        {Pid, Y} ->
            Y
    end.

aggregate_buffered_acc(_Pid, 0, Y) ->
    {continue, lists:reverse(Y)};
aggregate_buffered_acc(Pid, Count, Acc) ->
    receive
        {Pid, stop} ->
            {last, lists:reverse(Acc)};
        {Pid, Y} ->
            aggregate_buffered_acc(Pid, Count - 1, [Y | Acc])
    end.

aggregate_buffered(Pid, Count) ->
    aggregate_buffered_acc(Pid, Count, []).

calc_discrete_frequency(F, Fs) ->
    DiscreteFrequency = F / Fs,
    math:tan(?PI * DiscreteFrequency) / (2 * ?PI).

calc_discrete_cutting_frequency(Options) ->
    SamplingFrequency = proplists:get_value(sampling_frequency, Options),
    CuttingFrequency = proplists:get_value(cutting_frequency, Options),
    calc_discrete_frequency(CuttingFrequency, SamplingFrequency).

calc_discrete_ranged_frequency(Options) ->
    SamplingFrequency = proplists:get_value(sampling_frequency, Options),
    {F1, F2} = proplists:get_value(frequency_range, Options),
    Fd1 = calc_discrete_frequency(F1, SamplingFrequency),
    Fd2 = calc_discrete_frequency(F2, SamplingFrequency),
    {lists:min([Fd1, Fd2]), lists:max([Fd1, Fd2])}.

% Supported only two denominators.
gen_filter(lpf, Options) ->
    QualityFactor = proplists:get_value(cuality_factor, Options, 1 / math:sqrt(2)),
    Fc = calc_discrete_cutting_frequency(Options),
    Denominator = 1 + 2 * ?PI * Fc / QualityFactor
                  + 4 * math:pow(?PI, 2) * math:pow(Fc, 2),
    DelayerA = [(8 * math:pow(?PI, 2) * math:pow(Fc, 2) - 2) / Denominator,
                (1 - 2 * ?PI * Fc / QualityFactor
                 + 4 * math:pow(?PI, 2) * math:pow(Fc, 2)) / Denominator],
    DelayerB = [4 * math:pow(?PI, 2) * math:pow(Fc, 2) / Denominator,
                8 * math:pow(?PI, 2) * math:pow(Fc, 2) / Denominator,
                4 * math:pow(?PI, 2) * math:pow(Fc, 2) / Denominator],
    Pid = start_loop(DelayerA, DelayerB),
    Pid;
gen_filter(hpf, Options) ->
    QualityFactor = proplists:get_value(cuality_factor, Options, 1 / math:sqrt(2)),
    Fc = calc_discrete_cutting_frequency(Options),
    Denominator = 1 + 2 * ?PI * Fc / QualityFactor
                  + 4 * math:pow(?PI, 2) * math:pow(Fc, 2),
    DelayerA = [(8 * math:pow(?PI,2) * math:pow(Fc, 2) - 2) / Denominator,
                (1 - 2 * ?PI * Fc / QualityFactor
                 + 4 * math:pow(?PI, 2) * math:pow(Fc, 2)) / Denominator],
    DelayerB = [1 / Denominator,
                -2 / Denominator,
                1 / Denominator],
    Pid = start_loop(DelayerA, DelayerB),
    Pid;
gen_filter(bpf, Options) ->
    {Fl, Fh} = calc_discrete_ranged_frequency(Options),
    Denominator = 1 + 2 * ?PI * (Fh - Fl) + 4 * math:pow(?PI, 2) * Fl * Fh,
    DelayerA = [(8 * math:pow(?PI, 2) * Fl * Fh - 2) / Denominator,
                (1 - 2 * ?PI * (Fh - Fl) + 4 * math:pow(?PI, 2) * Fl * Fh) / Denominator],
    DelayerB = [2 * ?PI * (Fh - Fl) / Denominator,
                0,
                -2 * ?PI * (Fh - Fl) / Denominator],
    Pid = start_loop(DelayerA, DelayerB),
    Pid;
gen_filter(bef, Options) ->
    {Fl, Fh} = calc_discrete_ranged_frequency(Options),
    Denominator = 1 + 2 * ?PI * (Fh - Fl) + 4 * math:pow(?PI, 2) * Fl * Fh,
    DelayerA = [(8 * math:pow(?PI, 2) * Fl * Fh - 2) / Denominator,
                (1 - 2 * ?PI * (Fh - Fl) + 4 * math:pow(?PI, 2) * Fl * Fh) / Denominator],
    DelayerB = [(4 * math:pow(?PI, 2) * Fl * Fh + 1) / Denominator,
                (8 * math:pow(?PI, 2) * Fl * Fh - 2) / Denominator,
                (4 * math:pow(?PI, 2) * Fl * Fh + 1) / Denominator],
    Pid = start_loop(DelayerA, DelayerB),
    Pid;
gen_filter(resonance, Options) ->
    QualityFactor = proplists:get_value(cuality_factor, Options, 1 / math:sqrt(2)),
    Fc = calc_discrete_cutting_frequency(Options),
    Denominator = 1 + 2 * ?PI * Fc / QualityFactor
                  + 4 * math:pow(?PI, 2) * math:pow(Fc, 2),
    DelayerA = [(8 * math:pow(?PI, 2) * math:pow(Fc, 2) - 2) / Denominator,
                (1 - 2 * ?PI * Fc / QualityFactor
                 + 4 * math:pow(?PI, 2) * math:pow(Fc, 2)) / Denominator],
    DelayerB = [2 * ?PI * Fc / QualityFactor,
                0,
                -2 * ?PI * Fc / QualityFactor],
    Pid = start_loop(DelayerA, DelayerB),
    Pid;
gen_filter(notch, Options) ->
    QualityFactor = proplists:get_value(cuality_factor, Options, 1 / math:sqrt(2)),
    Fc = calc_discrete_cutting_frequency(Options),
    Denominator = 1 + 2 * ?PI * Fc / QualityFactor
                  + 4 * math:pow(?PI, 2) * math:pow(Fc, 2),
    DelayerA = [(8 * math:pow(?PI, 2) * math:pow(Fc, 2) -2) / Denominator,
                (1 - 2 * ?PI * Fc / QualityFactor
                 + 4 * math:pow(?PI, 2) * math:pow(Fc, 2)) / Denominator],
    DelayerB = [(4 * math:pow(?PI, 2) * math:pow(Fc, 2) + 1) / Denominator,
                (8 * math:pow(?PI, 2) * math:pow(Fc, 2) - 2) / Denominator,
                (4 * math:pow(?PI, 2) * math:pow(Fc, 2) + 1) / Denominator],
    Pid = start_loop(DelayerA, DelayerB),
    Pid.

%
% Followings are usage functions.
%

aggregate_filtered_acc(FilterPid, CallerPid, Acc) ->
    case aggregate_buffered(FilterPid, 512) of
        {last, Y} ->
            CallerPid ! {self(), lists:reverse(Y ++ Acc)};
        {continue, Y} ->
            aggregate_filtered_acc(FilterPid, CallerPid, Y ++ Acc)
    end.

aggregate_filtered(FilterPid, CallerPid) ->
    aggregate_filtered_acc(FilterPid, CallerPid, []).

apply_filter([], FilterPid, AggregaterPid) ->
    FilterPid ! {AggregaterPid, stop};
apply_filter([X | Xr], FilterPid, AggregaterPid) ->
    FilterPid ! {AggregaterPid, X},
    apply_filter(Xr, FilterPid, AggregaterPid).

filter(FilterName, Samples, Options) ->
    FilterPid = gen_filter(FilterName, Options),
    AggregaterPid = spawn_link(?MODULE, aggregate_filtered, [FilterPid, self()]),
    spawn_link(?MODULE, apply_filter, [Samples, FilterPid, AggregaterPid]),
    receive
        {AggregaterPid, Y} ->
            Y
    end.