信息处理 - 数字状态变量 (Chamberlin) 滤波器的频率响应 - 吾爱随笔录

数字状态变量 (Chamberlin) 滤波器的频率响应

信息处理无限脉冲响应频率响应

2022-02-23 13:10:42

我最近根据这里的建议实施了一个数字状态滤波器。我已经测试了该滤波器作为具有有限系数量化精度的非常低截止低通滤波器的性能，并且它完全按预期工作。但是，我只在时域中对其进行了测试。我还想根据它的频率响应来描述它。我使用Mason 的增益公式得出以下传递函数：

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

其中 $f=2\sin(\pi f_c/f_s)$ 和 $q=1/Q$ 是选择的常数（请参阅此链接）。根据该链接， $f_c=10\,\text{Hz}$ 是截止频率， $f_s=500\,\text{kHz}$ 是采样率， $Q=1/\sqrt{2}$ 。这是滤波器的框图（取自 Microprocessors 的音乐应用程序）以供参考：

但是，当我绘制响应（ $H(e^{j\omega})$ ）时，它看起来并不像我预期的那样。这是用于绘图的 Python 代码：

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()

这是绘制的频率响应

形状如我所料（低通和低 q 值，12dB/oct. 滚降）。但是，截止频率大约是远低于我设置的。此外，我对频率响应中的尖峰有些不安，这是我没想到的。我是否错误地设置了这个滤波器，或者错误地计算了频率响应？这是我第一次使用 Mason 的增益公式，所以我可能做错了。为什么我在的增益？如何获得正确的截止频率？这些“尖峰”是否值得关注？为什么它们会出现，我该如何移除它们？ $2\times 10^{-5}\,\text{Hz}$ $10\,\text{Hz}$ $-3\,\text{dB}$ $10\,\text{Hz}$

1个回答

您看到这些结果的主要原因是您的传递函数位于域中。独特的频率响应限于并且是周期性的。在您的情况下，您使用的无效值来定义域频谱存在的位置。这就是为什么你会得到这些尖峰。使用您的频率向量，我在下面复制了您的结果 $z$ $[-\pi, \pi]$ $f$ $z$

您生成的内容实际上是正确的，但是您以错误的方式看待它。

如果您选择的范围在之间，您将获得以下频率响应 $[0.1, \pi]$

它更简洁一些，您可以尝试如何定义函数和频率轴，使其看起来更理想。如果您查看数据标记，3-dB 点仍然是那个讨厌的 Hz。这仍然是正确的！只是这是一个“离散”频率。 $2 \times10^{-5}$

为了产生连续时间频率，您必须使用以下公式从离散域频率 $f$ $f'$

f = f^{'} f s Hz

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

使用我们得到的 3-dB 截止频率

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

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

这正是您所期望的。关于连续频率与离散频率以及 MATLAB 特定的详细信息，我有一个相关的答案。

编辑：使用过滤器

为了证明滤波器有效，我们生成并过滤信号

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

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

正弦曲线是、和。 $5 \text{ Hz}$ $10 \text{ Hz}$ $50 \text{ Hz}$

我使用了 MATLAB 的filter()函数，它采用传递函数的系数并生成一个差分方程来执行滤波。鉴于您的函数已经处于有理形式，识别系数是微不足道的。Python 应该有一个等价物。

下面是滤波前后 $x(t)$

您可以看到分量仍然存在，分量由于处于 3-dB 截止频率而被部分衰减，并且音调被抑制。下面是我用来生成这些结果的 MATLAB 代码。 $5 \text{ Hz}$ $10 \text{ Hz}$ $50 \text{ Hz}$

%% 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");

其它你可能感兴趣的问题

上一篇带差分方程和系统函数的单面 Z 变换下一篇基本时期x [ k n ]x[kn]