滤波器的画图、模拟信号与数字信号之间的关系、矩形窗滤波器的设计

本文主要讨论了数字滤波技术中各种符号的含义，以及模拟和数字信号之间的关系。

• 考虑位移变化
• 不考虑位移变化

从本质上说，一个时间采样的信号，所以采样率是fs；一个是以位移来采样的信号，如果采样间隔为L，那么采样率是1/L。

时域频率 Ω \Omega Ω与频率 f f f之间的关系是： Ω = 2 π f \Omega = 2\pi f Ω=2πf。

在模拟信号中，传输函数通常用于表示系统，因此传输函数中 s = j Ω = j 2 π f s = j\Omega = j2\pi f s=jΩ=j2πf，matlab验证如下：

傅立叶变换： F ( w ) = ∫ ? ∞ ∞ f ( t ) e ? i ω t d t F(w)=\int_{-\infty}^{ \infty} f(t) e^{-i \omega t} d t F(w)=∫ −∞+∞​ f(t)e−iωtdt 拉普拉斯变换： F ( ω ) = ∫ 0 + ∞ f ( t ) e − σ t e − i ω t d t = ∫ 0 + ∞ f ( t ) e − ( σ + i ω ) t d t F(\omega)=\int_{0}^{+\infty} f(t) e^{-\sigma t} e^{-i \omega t} d t=\int_{0}^{+\infty} f(t) e^{-(\sigma+i \omega) t} d t F(ω)=∫0+∞​f(t)e−σte−iωtdt=∫0+∞​f(t)e−(σ+iω)tdt

代码：

B = [1,2];
A = [3,4];
[ H , tmp ] = freqs(B, A);%# 连续的模拟信号
figure;
semilogx(tmp/pi/2 , 20*log10(H),'--','LineWidth',1);
set(gca,'Fontname','Times New Roman','fontsize',16);
xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
f = 0.1:0.001:100;
Omega = 2*pi*f;
Cs = (1j*Omega*B(1) + B(2) )./( 1j*Omega*A(1) + A(2) );
hold on
semilogx( f , 20*log10(Cs),'LineWidth',1 );
xlabel('模拟频率 (Hz)');
ylabel('dB');   

总之这里注意2*pi的坐标轴变化，具体原因不了解。

在数字信号中，常常用 ω \omega ω表示数字频率， T i T_i Ti​表示相邻采样点之间的时间间隔，则有： ω = Ω T i = 2 π f / f s \omega = \Omega T_i = 2\pi f/f_s ω=ΩTi​=2πf/fs​，另外常常用 z z z表示数字系统特性， z = e x p ( i ω ) = e x p ( i 2 π f / f s ) z = exp(i\omega) = exp(i2\pi f/f_s) z=exp(iω)=exp(i2πf/fs​)。

psi = 0:0.001:900;
fs = 1000;
wd = 0.001;
b = [1-wd,-1+wd];
a = [1,-1+wd];
z_1 = exp(-j*2*pi*psi/fs);
Rz = (1-wd)*(1 - z_1)./( 1-(1-wd) * z_1 );
figure;semilogx(psi,20*log10(Rz),'LineWidth',1);
[h,f] = freqz(b,a,1e5,fs);
hold on;semilogx(f,20*log10(h),'--r','LineWidth',2)

一般模拟滤波器的传递函数为 H a ( s ) H_a(s) Ha​(s)，其时域频率特性为 H z ( j Ω ) H_z(j\Omega) Hz​(jΩ)。数字滤波器的转移函数为 H d ( z ) H_d(z) Hd​(z)，其频率特性可以表示为 H d ( e i ω ) H_d(e^{i\omega}) Hd​(eiω)。

时间域频率 Ω \Omega Ω与空间域频率 ψ \psi ψ之间的关系为： Ω = 2 π ψ v \Omega = 2\pi \psi v Ω=2πψv。

数字频率： ω = Ω T i = 2 π ψ ( v T i ) = 2 π ψ 0.25 \omega = \Omega T_i = 2\pi \psi (vT_i) = 2\pi \psi 0.25 ω=ΩTi​=2πψ(vTi​)=2πψ0.25。（这里假定0.25m采样）

所以： z = e x p ( 2 π ψ 0.25 ) z = exp(2\pi \psi 0.25) z=exp(2πψ0.25)。

那么这里需要注意的是：时间频率 f f f与空间频率 ψ \psi ψ的关系是： f = v ψ f = v\psi f=vψ。同时，空间域频率 ψ = 1 / λ \psi = 1/\lambda ψ=1/λ，这个条件很重要。

因为是基于空间域采样，每隔0.25m采样一次，那么相当于采样率为4Hz。

矩形窗是设计滤波器的常用方法，下面主要考虑基于位移移变得情况下的滤波器设计。

矩形窗窗函数可表示为： w [ k ] = { 1 − M ≤ k ≤ M 0  other  w[k]=\left\{\begin{array}{cc} 1 & -M \leq k \leq M \\ 0 & \text { other } \end{array}\right. w[k]={ 10​−M≤k≤M other ​ 矩形窗滤波器系统函数可表示为： W ( z ) = 1 2 M + 1 z M − z − M − 1 1 − z − 1 W(z)=\frac{1}{2 M+1} \frac{z^{M}-z^{-M-1}}{1-z^{-1}} W(z)=2M+11​1−z−1zM−z−M−1​ 幅度函数： W ( Ω ) = 1 N sin ⁡ ( N Ω / 2 ) sin ⁡ ( Ω / 2 ) , N = 2 M + 1 W(\Omega)=\frac{1}{N} \frac{\sin (N \Omega / 2)}{\sin (\Omega / 2)}, \quad N=2 M+1 W(Ω)=N1​sin(Ω/2)sin(NΩ/2)​,N=2M+1 以下的代码基于M=30的条件进行仿真。

• 幅度函数

M = 30;
N = 2*M+1;
lamda = 1:0.1:1000;
pesi  = 1./lamda;           %空间域频率
omega = 2*pi*pesi*0.25;
Wn = sin(N.*omega/2)./sin(omega./2)./N;
hold on;semilogx(lamda,20*log10(abs(Wn)));

• 系统函数

K = 30;
M = 2*K+1;
lamda = 1:0.1:1000;
    pesi  = 1./lamda;           %空间域频率
    temp = exp(1j*2*pi.*pesi*0.25);%这个很重要
    %%换一种表示，就是两个矩形窗
    DjpesiDz =     ( temp.^K - temp.^(-K-1) )...
    ./( ( 1- temp.^(-1) ) ) ...
    /M ;
    hold on;semilogx( lamda , 20*log10(DjpesiDz) );

• freqz函数

b = ones(1,N)/N;
[h,f] = freqz(b,[zeros(1,N-1),1],10000,4);
hold on;semilogx(1./f,20*log10(h));grid on;

可见矩形窗滤波器对于信号相位的影响还是很大的，并且经常会产生信号相位的跳变。

理论上来说矩形窗应该是线性相位的。

以积分函数1/s举例，用反向差分法将其转换成为z传递函数。就可以看出来用数字系统去逼近模拟系统的时候会带来的问题，就是在某些频率的幅值上会做出改变。

幅值和相位在高频的部分都发生了改变。这个主要的原因就是本质上反向差分法是利用一阶泰勒展开逼近积分函数，这一定会忽略掉高阶项导致精度损失。高阶项对应的便可以认为是高频部分。

显然T的取值影响到最后的取值结果，假设采样率为4Hz（与轨检相同），那么有：

竖线的位置正好是4Hz，这与采样定理有关系。

这能说明什么？对于积分来说，对于越是高频（短波长）的处理，会导致更大的误差？可以这样理解？

双线性变换就是梯形积分？黑色是双线性变换的结果，只能说两者不一样，梯形积分可能会带来更多的误差。