# ModalRes

This is a `sndkit`

algorithm. A more up-to-date version
can be found here.

### 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.

*<<modalres.c>>=*```
#include <math.h>
#define SK_MODALRES_PRIV
#include "modalres.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
<<funcs>>
```

*<<modalres.h>>=*```
#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.

*<<funcdefs>>=*`void sk_modalres_init(sk_modalres *mr, int sr);`

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

### Struct

The encapsulating struct is called `sk_modalres`

.

*<<typedefs>>=*`typedef struct sk_modalres sk_modalres;`

*<<structs>>=*```
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.

*<<sk_modalres>>=*```
SKFLT x;
SKFLT y[2];
```

*<<init>>=*```
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.

*<<sk_modalres>>=*```
SKFLT b1;
SKFLT a1;
SKFLT a2;
```

*<<init>>=*```
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.

*<<sk_modalres>>=*`int sr;`

*<<init>>=*`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`

.

*<<sk_modalres>>=*`SKFLT s;`

*<<init>>=*`mr->s = 0;`

### Parameters

#### Frequency

Set with `sk_modalres_freq`

.

*<<funcdefs>>=*`void sk_modalres_freq(sk_modalres *mr, SKFLT freq);`

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

*<<sk_modalres>>=*```
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.

*<<init>>=*```
sk_modalres_freq(mr, 1000);
mr->pfreq = -1;
```

#### Q value

Set with `sk_modalres_q`

.

*<<funcdefs>>=*`void sk_modalres_q(sk_modalres *mr, SKFLT q);`

*<<funcs>>=*```
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.

*<<sk_modalres>>=*```
SKFLT q;
SKFLT pq;
```

*<<init>>=*```
sk_modalres_q(mr, 1);
mr->pq = -1;
```

### Computation

A single sample is computed with `sk_modalres_tick`

.

*<<funcdefs>>=*`SKFLT sk_modalres_tick(sk_modalres *mr, SKFLT in);`

*<<funcs>>=*```
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}}
```

*<<update_coefficients>>=*```
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.

*<<compute_difference_equation>>=*`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`

.

*<<update_filter_state>>=*```
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.

*<<scale_output>>=*`out *= mr->s;`