ModalRes

Overview

ModalRes is an implementation of a 2nd-order resonator filter suitable for use in modal synthesis.

This filter is particularly useful in constructing bell or other impact-based sounds. In order to do this, you'd need to know the modal ratios pre-measured of time. This is what is known as black-box modelling. A similar example of this is formant-based speech synthesis.

Some good modal ratios and measurements can be found here.

Tangled Files

modalres.c and modalres.h are the two files that are tangled.

#include <math.h>
#define SK_MODALRES_PRIV
#include "modalres.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

<<funcs>>

#ifndef SK_MODALRES_H
#define SK_MODALRES_H

#ifndef SKFLT
#define SKFLT float
#endif

<<typedefs>>
<<funcdefs>>

#ifdef SK_MODALRES_PRIV
<<structs>>
#endif
#endif

Initiazation

Initialized with sk_modalres_init. Only the sampling rate needs to be supplied.

void sk_modalres_init(sk_modalres *mr, int sr);

void sk_modalres_init(sk_modalres *mr, int sr)
{
    <<init>>
}

Struct

The encapsulating struct is called sk_modalres.

typedef struct sk_modalres sk_modalres;

struct sk_modalres {
    <<sk_modalres>>
};

Filter State Variables and Coefficients

The filter state parameters are the typical things you'd need in a biquadratic filter difference equation.

Variables are created to store the previous x value and 2 y values.

SKFLT x;
SKFLT y[2];

mr->x = 0;
mr->y[0] = 0;
mr->y[1] = 0;

Coefficients for the biquad are b1 for the x value, and a1 and a2 for the y values. These coefficents get computed in the computation section.

SKFLT b1;
SKFLT a1;
SKFLT a2;

mr->b1 = 0;
mr->a1 = 0;
mr->a2 = 0;

As per usual, a copy of the sampling rate is stored, as it is needed to calculate filter coefficients.

int sr;

mr->sr = sr;

An unusual part of this filter design is the scaler variable, which is used to scale the audio output, but is not part of the actual transer function. We'll call this one s.

SKFLT s;

mr->s = 0;

Parameters

Frequency

Set with sk_modalres_freq.

void sk_modalres_freq(sk_modalres *mr, SKFLT freq);

void sk_modalres_freq(sk_modalres *mr, SKFLT freq)
{
    mr->freq = freq;
}

SKFLT freq;
SKFLT pfreq;

The frequency parameter uses caching to check if filter coefficients needs to be updated. The cached value is set to be negative in order to force computation on the first tick.

sk_modalres_freq(mr, 1000);
mr->pfreq = -1;

Q value

Set with sk_modalres_q.

void sk_modalres_q(sk_modalres *mr, SKFLT q);

void sk_modalres_q(sk_modalres *mr, SKFLT q)
{
    mr->q = q;
}

The Q parameter uses caching to check if filter coefficients needs to be updated. The cached value is set to be negative in order to force computation on the initial tick.

SKFLT q;
SKFLT pq;

sk_modalres_q(mr, 1);
mr->pq = -1;

Computation

A single sample is computed with sk_modalres_tick.

SKFLT sk_modalres_tick(sk_modalres *mr, SKFLT in);

SKFLT sk_modalres_tick(sk_modalres *mr, SKFLT in)
{
    SKFLT out;

    out = 0;
    <<update_coefficients>>
    <<compute_difference_equation>>
    <<update_filter_state>>
    <<scale_output>>
    return out;
}

Update the coefficients, if needed. First, frequency is converted to radians. While traditionally represented as the lowercase $\omega$ in mathematical notation, the variable w will be used in C as crude approximation.

From there the following coefficients are computed:

b1 = {1 \over \omega^2 + {\omega \over 2Q}}
a1 = {1 - 2\omega^2}b1
a2 = {\omega^2 - {\omega \over 2Q}}

if (mr->freq != mr->pfreq || mr->q != mr->pq) {
    SKFLT w;
    SKFLT a, b, d;

    w = mr->freq * 2.0 * M_PI;

    a = mr->sr / w;
    b = a*a;
    d = 0.5*a;

    mr->pfreq = mr->freq;
    mr->pq = mr->q;

    mr->b1 = 1.0 / (b + d/mr->q);
    mr->a1 = (1.0 - 2.0*b) * mr->b1;
    mr->a2 = (b - d/mr->q) * mr->b1;
    mr->s = d;
}

Compute difference equation. This computes the difference equation:

$$
y(n) = b0 x(n) + b1 x(n - 1) - a1 y(n - 1) - a2 y(n - 2)
$$

Where $b0$ is 0, effectively cancelling out this term.

What's interesting about this is that it adds an implicit one-sample delay to the filter. This is done because more accurately simulates the impulse response in an ideal mass-spring-damper system. The produces is, after all, a response to the impulse.

out = mr->b1*mr->x - mr->a1*mr->y[0] - mr->a2*mr->y[1];

Update filter state. Shift things forward in time one sample. What was once y[0] is now y[1]. What is now the current output becomes y[0], what is now the current input value because the previous input value x.

mr->y[1] = mr->y[0];
mr->y[0] = out;
mr->x = in;

The output is scaled by the variable s before it is returned. It is not part of the recursive filter function.

out *= mr->s;