Frequency response of a digital state variable (Chamberlin) filter

Question

I recently implemented a digital state filter based on the recommendation here. I've tested this filter's performance as a very low cutoff low-pass filter with limited coefficient quantization precision and it's working entirely as expected. However, I've only tested it in the time domain. I'd additionally like to characterize it in terms of its frequency response. I've used Mason's gain formula to arrive at the following transfer function:

$ H(z) = \frac{f^2z^{-1}}{1-z^{-1}(2-qf-f^2)+z^{-2}(1-qf)}, $

where $f=2\sin(\pi f_c/f_s)$ and $q=1/Q$ are chosen constants (see this link). According to that link, $f_c=10\,\text{Hz}$ is the cutoff frequency, $f_s=500\,\text{kHz}$ is the sampling rate, and $Q=1/\sqrt{2}$. Here's a block diagram of the filter (taken from Musical Applications of Microprocessors) for reference:

However, when I plot the response ($H(e^{j\omega})$) it doesn't quite look as I expect. Here's the Python code for plotting:

import numpy as np
import matplotlib.pyplot as plt

fc = 1e1
fsample = 500e3
fnyquist = fsample / 2
q = 1 / np.sqrt(2)
Fc = 2 * np.sin(np.pi * fc / fsample)
Q = 1 / q

def tf(f):
    w = 2 * np.pi * f
    z = np.exp(-1 * 1j * w)
    return (
        Fc ** 2
        * z
        / (1 - z * (2 - Q * Fc - Fc ** 2) + z ** 2 * (1 - Q * Fc))
    )

freq = np.logspace(-10, np.log10(fnyquist), int(1e5))
resp = [20 * np.log10(abs(tf(f))) for f in freq]
_, ax = plt.subplots()
ax.plot(freq, resp)
ax.grid(b=True, which="major")
ax.set_ylim(-120, 10)
ax.set_xscale("log")
plt.show()

Here's the plotted frequency response

The shape is as I would have expected (low-pass and low q-value with 12dB/oct. rolloff). However, the cutoff frequency, which is roughly $2\times 10^{-5}\,\text{Hz}$ is much lower than the $10\,\text{Hz}$ I set. Additionally, I'm somewhat perturbed by the spikes in the frequency response, which I didn't expect. Have I set up this filter incorrectly, or calculated the frequency response incorrectly? This is my first time using Mason's gain formula, so it's possible I've done that incorrectly. Why do I not see a gain of $-3\,\text{dB}$ at $10\,\text{Hz}$? How can I achieve the correct cutoff frequency? Are those "spikes" a cause for concern? Why are they present and how can I remove them?

Answer 1

The main reason you're seeing these results is because your transfer function is in the $z$-domain. The unique frequency response is limited to $[-\pi, \pi]$ and is periodic. In your case, you are using invalid values of $f$ to define where the $z$-domain spectrum exists. This is why you get those spikes. Using your frequency vector, I've replicated your result below

What you're generating is somewhat actually correct, but you're looking at it in the wrong way.

If you choose the range to be between $[0.1, \pi]$, you'll get the following frequency response

It's a little cleaner and you can play around with how you define your function and frequency axis to make it look more ideal. If you look at the data marker, the 3-dB point is still at that pesky $2 \times10^{-5}$ Hz. This is still correct! It's just that this is a "discrete" frequency.

In order to yield the continuous-time frequency $f$, you must convert from the discrete-domain frequency $f'$ by using the formula

$$f = f'fs \space \text{Hz}$$

Using the 3-dB cutoff frequency we get

$$f = (2 \times 10^{-5})(500 \times 10^3 )= 10 \space \text{Hz}$$

Which is exactly what you expect. There's a related answer I have regarding continuous vs discrete frequencies along with MATLAB-specific details here.

EDIT: Using the filter

To show that the filter works, we generate and filter the signal

$$x(t) = \cos(2\pi(5)t) + \cos(2\pi(10)t) + \cos(2\pi(50)t)$$

Which are sinusoids are $5 \text{ Hz}$, $10 \text{ Hz}$, and $50 \text{ Hz}$.

I used MATLAB's filter() function which takes the transfer function's coefficients and generates a difference equation to perform the filtering. Identifying the coefficients is trivial given that your function is already in rational form. Python should have an equivalent.

Below is the signal $x(t)$ before and after filtering

You can see that the $5 \text{ Hz}$ component survives, the $10 \text{ Hz}$ component is partially attenuated since it's at the 3-dB cutoff, and the $50 \text{ Hz}$ tone is suppressed. Below is the MATLAB code I used to generate these results.

%% Sampling and constants

fc = 10;
fs = 500e3;
fn = fs/2;

q = 1/sqrt(2);
Q = 1/q;
Fc = 2*sin(pi*fc/fs);

%% Manually define the transfer function. Uncomment to generate and manually plot the frequency response.
% f = logspace(-10, pi, 1e5);
% w = 2.*pi.*f;
% z = exp(-1i.*w);
% 
% freqResponse = (Fc.^2.*z)./(1 - z.*(2 - Q.*Fc - Fc.^2) + z.^2.*(1 - Q.*Fc));
% 
% figure;
% semilogx(f, 20*log10(abs(freqResponse)));
% xlabel("Normalized Frequency (Hz/sample)");
% ylabel("Magnitude (dB)");
% axis tight;
% ylim([-120 10]);

%% Using built-in function filter()

b = [0 Fc.^2];
a = [1 -(2 - Q.*Fc - Fc.^2) (1 - Q.*Fc)];
[h, w] = freqz(b, a, 1e5);

figure;
semilogx(w./(2*pi), 20*log10(abs(h)));
axis tight;
ylim([-120 10]);

%% Use the filter to process a signal

t = 0:1/fs:2;

x = cos(2*pi*(5).*t) + cos(2*pi*(10).*t) + cos(2*pi*(50).*t);

nfft = 10*numel(x);
f = fs.*(-nfft/2:nfft/2-1)./nfft;

figure;
subplot(2, 1, 1);
plot(f, abs(fftshift(fft(x, nfft)./nfft)).^2);
xlim([-100 100])
xlabel("Frequency (Hz)");
ylabel("Magnitude");
title("Original Signal");

subplot(2, 1, 2);
plot(f, abs(fftshift(fft(filter(b, a, x), nfft)./nfft)).^2);
xlim([-100 100])
xlabel("Frequency (Hz)");
ylabel("Magnitude");
title("Filtered Signal");