滤波器时间域与空间域幅频特性matlab绘图方法

1. 模拟滤波器

首先，明确模拟时间（角）频率与空间频率的关系：

$$\omega =2\pi v/\lambda =2\pi \varphi v$$

式中：

$\omega$为模拟频率，rad/s

$v$为速度，m/s

$\lambda$为波长，m

$\varphi$为空间频率，1/m

以二阶巴特沃斯滤波器为例。

系统函数为：

$$H(s)=\frac{a^2}{s^2+as+a^2}$$

使用matlab程序画出滤波器幅频特性，代码如下：

clear;
% H(s) = a^2 / (s^2 + a*s + a^2)
a = 10;
n = a^2;
d = [1 a a^2];
v = 10;

omega = logspace(-5,5,10000);
phi = omega ./ (2*pi*v);
lambda = 1./phi;

mag = freqs(n,d,omega);
mag = 20*log10(abs(mag));

figure;
semilogx(omega,mag,'LineWidth',2);hold on;grid on;
xlabel('1/s');ylabel('dB');title('滤波器时间域幅频特性');
axis([10e-2 10e2 -80 20]);
figure;
semilogx(phi,mag,'LineWidth',2);hold on;grid on;
xlabel('1/m');ylabel('dB');title('滤波器空间域幅频特性1');
axis([10e-3 10e1 -80 20]);
figure;
semilogx(lambda,mag,'LineWidth',2);hold on;grid on;
xlabel('m');ylabel('dB');title('滤波器空间域幅频特性2');
axis([10e-2 10e2 -80 20]);

上述程序中，思路是使用freqs()函数求出相应的频率响应。画图时(x,y)相同的是y，不同的是x，替换相应的x即可。本例程使用的参数是：

a = 10rad/s

v = 10m/s

2. 数字滤波器

同样，需要明确以下关系：

$$\Omega =\omega T=2\pi \Delta x/\lambda =2\pi \varphi \Delta x$$

式中：

$\Omega$为数字频率，rad

$\Delta x$为空间采样频率，m

以滑动平均滤波器为例。

系统函数为：

$$H(z)=\frac{1}{2M+1}\frac{z^M-z^{-M-1}}{1-z^{-1}}$$

可以求得其频谱函数为：

$$W(e^{j\Omega })=\frac{sin(N\Omega /2)}{Nsin(\Omega /2)},N=2M+1$$

代码如下：

clear;

M = 10;
N = 2*M+1;
a = ones(1,N)/N;

deltax = 0.1;
Omega = logspace(-5,5,10000);
phi = Omega ./ (2*pi*deltax);
lambda = 1./phi;

mag_freqz = freqz(a,1,Omega);
mag_freqz = 20*log10(abs(mag_freqz));

mag_direct = sin(N*Omega/2)./sin(Omega/2) / N;
mag_direct = 20*log10(abs(mag_direct));

figure;
semilogx(Omega,mag_freqz,'LineWidth',2);hold on;grid on;
semilogx(Omega,mag_direct-10,'LineWidth',2);hold on;grid on;
xlabel('rad');ylabel('dB');title('滤波器时间域幅频特性');
axis([10e-5 10e0/4 -80 20]);
figure;
semilogx(phi,mag_freqz,'LineWidth',2);hold on;grid on;
semilogx(phi,mag_direct-10,'LineWidth',2);hold on;grid on;
xlabel('1/m');ylabel('dB');title('滤波器空间域幅频特性1');
axis([10e-5 10e0/2 -80 20]);
figure;
semilogx(lambda,mag_freqz,'LineWidth',2);hold on;grid on;
semilogx(lambda,mag_direct-10,'LineWidth',2);hold on;grid on;
xlabel('m');ylabel('dB');title('滤波器空间域幅频特性2');
axis([10e-2*2 10e3 -80 20]);

上述程序中，使用离散序列ones(1,N)/N代替了系统函数的表达式。

分别使用freqz()方法和频谱函数方法求其幅频特性。

绘图方法与上面方法相同，在20、25、30行中人为设置了曲线的偏移，去掉偏移后，两种方法绘图完全一致。