现代雷达系统分析与设计---动目标显示(MTI)

MTI是指利用杂波抑制滤波器来抑制各种杂波，提高雷达信号的信杂比，以利于运动目标检测的技术。以地杂波为例，杂波谱通常集中于直流(多普勒频率$f_d=0$)和雷达重复频率$f_r$的整数倍处，如下图(a)所示。在连续波雷达中，由于多数情况下杂波功率集中于零频附近，因此杂波可以通过忽略直流输出来避免或抑制。在脉冲雷达中MTI滤波器就是利用杂波与运动目标的多普勒频率的差异，使得滤波器的频率响应在杂波谱的位置形成"凹口"，以抑制杂波，而让动目标回波通过后的损失尽量小或没有损失。为了有效地抑制杂波，MTI滤波器需要在直流和PRF的整数倍处具有较深的阻带。下图(b)显示了一个典型的MTI滤波器的频率响应图，下图©显示了输入下图(a)所示的功率谱密度时的滤波器输出。

延迟线对消器

延迟线对消器是最早出现，也是最常用的MTI滤波器之一。根据对消次数的不同，又分为单延迟线对消器、双延迟线对消器和多延迟线对消器。

单延迟线对消器

单延迟线对消器如下图所示。它由延迟时间等于发射脉冲重复周期PRI($T_r$)的延迟单元(数字延迟线)和加法器组成。单延迟线对消器经常称为"两脉冲对消器"或者"一次对消器"。

对消器的脉冲响应表示为$h(t)$，输出$y(t)$等于脉冲响应$h(t)$与输入$x(t)$之间的卷积。

输出信号$y(t)$为

$y(t)=x(t)-x(t-T_r)$

对消器的脉冲响应为

$h(t)=\delta (t)-\delta (t-T_r)$

式中，$\delta (\bullet )$为$\delta$函数。由此可以得到$h(t)$的傅里叶变换(FT)，即频率响应为

$H(\omega )=1-e^{-j\omega T_r}$

式中，$\omega =2\pi f$

在z域，单延迟线对消器的传递函数为

$H(z)=1-z^{-1}$

单延迟线对消器的功率增益为

$|H(\omega ){|}^2=H(\omega )H^{\ast }(\omega )=(1-e^{-j\omega T_r})(1-e^{j\omega T_r})=2(1-cos\omega T_r)=4(sin(\frac{\omega T_r}{2}){)}^2$

双延迟线对消器

双延迟线对消器如下图所示。它由两个单延迟线对消器级联而成。双延迟线对消器经常称为"三脉冲对消器"或者"二次对消器"。

双延迟线对消器脉冲响应为

$h(t)=\delta (t)-2\delta (t-T_r)+\delta (t-2T_r)$

双延迟线对消器的功率增益为

$|H(\omega ){|}^2=|H_1(\omega ){|}^2|H_1(\omega ){|}^2$

式中，$|H_1(\omega ){|}^2$为单延迟线对消器的功率增益。由此得到

$|H(\omega ){|}^2=16(sin(\frac{\omega T_r}{2}){)}^4$

在z域，其传递函数为

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

下图给出了单延迟线对消器和双延迟线对消器的归一化频率响应。

从图中可以看出，双延迟线对消器比单延迟线对消器具有更好的响应(更深的凹口和更平坦的通带响应)。单延迟线对消器的频率响应较差，原因在于其阻带没有宽的凹口。而双延迟线对消器无论在阻带还是在通带上都比单延迟线对消器有更好的频率响应，因此比单延迟线对消器得到更广泛的应用。

单延迟线对消器是一个非常简单的滤波器。它的实现不需要乘法运算，每一个输出采样只需要一次减法运算。然而，与理想的高通滤波器相比，它是一个很差的近似。双延迟线对消器能够明显地改善零多普勒附近地凹口宽度，但并不能改善非零多普勒频率处地频率响应。双延迟线对消器的每一个输出采样只需要两次减法运算。

多延迟线对消器

多延迟线对消器是由多个单延迟线对消器级联而成，N延迟线对消器的脉冲响应为

$h(t)={\sum }_{n=0}^N{\omega }_n\delta (t-n{T}_r)$

式中，N为对消器的次数，对消器的系数${\omega }_n$为二项式系数，用下式计算

${\omega }_n=(-1{)}^n{C}_N^n=(-1{)}^n\frac{N!}{(N-n)!n!}，n=0,1,\dots \dots ,N$

所以，多延迟线对消器可用下图来统一表示

N次对消器传递函数为

$H(z)=(1-z^{-1}{)}^N={\sum }_{n=0}^N{\omega }_n{z}^{-n}$

从上式可以看出，N次对消器在$z=1$处有N重零点。N次对消器的频率响应为

$H(\omega )=(1-e^{-j\omega T_r}{)}^N$

其幅频响应和相频响应分别为

$|H(\omega )|=|2sin(\frac{\omega T_r}{2}){|}^N$

$\varphi (\omega )=N(\frac{\pi }{2}-\frac{\omega T_r}{2})$

可见，相位响应$\varphi (\omega )$与$\omega$是线性关系。所以对消器是一种线性相位滤波器，回波信号通过它后，相位关系不产生非线性变化。

假设输入杂波具有中心频率为零的高斯型功率谱，其功率谱的标准偏差为${\delta }_f$，则对消器的改善因子仅与${\delta }_f$和雷达重复频率$f_r$有关。N脉冲MTI改善因子的通用表达式为

$I_N=\frac{{\sum }_{n=0}^{N-1}{\omega }_n^2}{(2(N-1)-1)!!}(\frac{f_r}{2\pi {\delta }_f}{)}^{2(N-1)}=\frac{Q^2}{(2(N-1)-1)!!}{K}^{-2(N-1)}$

式中，$Q^2={\sum }_{n=0}^{N-1}{\omega }_n^2$，为MTI滤波器的二项式系数的平方和；$K=2\pi {\delta }_f/f_r$；双阶乘符号的定义为

$(2N-1)!!=1\ast 3\ast 5\ast \dots \dots \ast (2N-1)，0!!=1$

N脉冲MTI对消器的系数${\omega }_n$和改善因子如下表。

下图给出了改善因子与归一化谱宽$({\delta }_f/f_r)$之间的关系。可见，改善因子主要取决于杂波谱的归一化谱宽。

MTI滤波器如果具有与杂波功率谱主峰宽度相适应的滤波凹口和相对平坦的通带，会使杂波抑制滤波器输出的目标信号在通带内不随$f_d$而变化。但是延迟线对消器不能满足这种要求，当目标的多普勒频率为雷达重复频率的整数倍时，目标将被对消掉。为了解决这个问题，需要采用脉间重复频率参差的方法，并优化设计MTI滤波器。

参差重复频率

所有离散时间滤波器的频率响应都是周期性的，其周期在归一化频率域为单位1。既然MTI滤波器设计在零多普勒频率处形成凹口，那它们也同样会在整数倍$f_r$的多普勒频率处形成凹口。因此，当运动目标的多普勒频率等于整数倍$f_r$时，这些运动的目标也会被MTI滤波器滤掉。对应于这些多普勒频率(即脉冲重复频率的整数倍)的径向速度称为盲速，这是因为具有这些径向速度的运动目标也会被抑制掉。也就是说，系统对于这些目标是"盲"的。从数字信号处理的角度看，盲速代表那些模糊到零多普勒频率的目标速度。参差重复频率是一种可以用来防止盲速影响的措施。

盲速

对于发射脉冲重复频率为$f_r$的脉冲雷达，如果运动目标相对雷达的径向速度$v_r$引起的相邻周期回波信号相位差$\Delta \varphi =2\pi f_d{T}_r$，其中，$f_d=2v_r/\lambda$为$v_r$产生的多普勒频率，$T_r$为雷达脉冲重复周期。当$\Delta \varphi$为$2\pi$的整数倍时，由于脉冲雷达系统对目标多普勒频率取样的结果，相位检波器的输出为等幅脉冲，与固定目标相同，因此动目标显示输出为零，这时的目标速度称为盲速，具体推导如下

$\Delta \varphi =2\pi f_{bn}{T}_r=n2\pi$

式中，$f_{bn}$为产生盲速时的目标多普勒频率。

$f_{bn}=\frac{n}{T_r}=n{f}_r$

所以，盲速$v_{bn}$为

$v_{bn}=\frac{1}{2}n\lambda f_r$

对于一个给定的$f_r$，不模糊距离为$R_u=c/(2f_r)$。当增加$f_r$时，距离的不模糊范围$R_u$减小，而第一个盲速增加。

盲速可以通过选择足够高的PRF来避免，因为当PRF足够高时，可以使第一个盲速超过任何可能的实际目标速度。然而，遗憾的是，较高的PRF对应于较短的不模糊距离。通常情况下，无法得到一个PRF能够同时满足要求的不模糊距离和多普勒覆盖区。

为了解决盲速的问题，常用的方法是利用参差重复频率，它能大大提高第一个盲速，而不会减小不模糊距离。PRF参差的实现既可以是脉间参差，也可以脉组参差。

脉间参差的优点是能够在一个驻留时间内提高不模糊的多普勒覆盖区。但脉间参差的一个缺点是数据是非均匀采样的序列，这使得应用相干多普勒处理变得困难，而且也使分析变得复杂。另一个缺点是模糊的主瓣杂波会导致脉冲间得杂波幅度随着PRF的变化而变化，这是由于距离模糊的杂波(来自前面的脉冲)会随着PRF的变化而折叠到不同的距离单元中去。因此，脉间PRF参差通常只用于无距离模糊的低PRF工作模式。

参差重复频率

如果雷达采用N个重复频率$f_{r1}，f_{r2}，\dots ，f_{rN}$，它们的重复周期可以表示为

$$\{\begin{array}{c}{T}_{r1}=1/f_{r1}=K_1\Delta T\\ T_{r2}=1/f_{r2}=K_2\Delta T\\ \dots \dots \\ T_{rN}=1/f_{rN}=K_N\Delta T\end{array}$$

式中，$\Delta T$为$[T_{r1},T_{r2},\dots \dots ，T_{rN}]$的最大公约周期，则参差周期之比为

$T_{r1}：T_{r2}：\dots \dots ：T_{rN}=K_1：K_2：\dots \dots ：K_N$

式中，$[K_1：K_2：\dots \dots ：K_N]$为参差码，参差码中最大K值与最小K值之比称为参差周期的最大变比$r$，

$r=max[K_1：K_2：\dots \dots ：K_N]/min[K_1：K_2：\dots \dots ：K_N]$

如果$K_i$之间互异互素，则第一个真正的盲速对应的多普勒频率$f_{bn}$为

$f_{bn}=\frac{1}{\Delta T}$

雷达的平均重复周期为

$T_r=\frac{1}{N}{\sum }_{i=1}^N{T}_{ri}=K_{av}\Delta T$

式中，$K_{av}$为参差码的均值。因此

$K_{av}=\frac{T_r}{\Delta T}=T_r{f}_{bn}=\frac{f_{bn}}{f_r}$

$f_{bn}=K_{av}{f}_r$

因为$f_r=1/T_r$是雷达平均重复频率，所以也称$K_{av}$为盲速扩展倍数

参数PRF时的三脉冲对消器如下图所示。

每个脉冲之间MTI滤波器的系数不同，因此，它是一种时变滤波器。如果雷达依次采用三种重复频率$T_1、T_2、T_3、T_1，\dots \dots$(常称为"三变T"),则有三组MTI滤波器依次工作。

参差MTI滤波器的频率响应取决于参差周期和滤波器权矢量。如果滤波器权矢量为二项式系数，就构成了参差对消器。

代码：

len = 3; %滤波器长度
bianT = [27 28 29]; %变T码
f = (-1:0.01:28); %滤波器归一化响应频率范围
fr = 1/mean(bianT);
bianT_num = length(bianT);
bianT = repmat(bianT,1,len);
T = zeros(1,len);
f_ran = f .* fr; %频率范围
hd = zeros(bianT_num, length(f)); %频率响应
N = len - 1;
for m = 0:N
    w(m+1) = (-1)^(m-1)*factorial(N)/(factorial(m)*factorial(N-m)); %计算权系数
end
w = w./max(abs(w)); %归一化权系数
for k = 1:bianT_num
    for m = 2:len
        T(m) = sum(bianT(k:k+m-2));
    end
    hd(k,:) = w*exp(-1i*2*pi*T'*f_ran);
end
Hd = 20*log10(abs(hd));
figure;
plot(f_ran./fr,Hd);

运行结果：

从上图中可以看出，使用参差重复频率能在很大程度上提高第一盲速

参数MTI滤波器速度响应凹口的深度与对消器的形式无关，也与雷达天线波束内接收的脉冲数无关，而只和参差周期的最大变比有关。最大变比越大，对应的凹口深度越浅。

优化MTI滤波器

滤波器主要分为无限脉冲响应(IIR)滤波器和有限脉冲响应(FIR)滤波器。IIR滤波器的优点是可用相对较少的阶数达到预期的滤波器响应，但是其相位特性是非线性的，在MTI滤波器中很少采用。而FIR滤波器具有线性相位特性，所以MTI滤波器主要采用FIR滤波器。延迟线对消器也是一种FIR滤波器，是系数符合二项式展开式的特殊FIR滤波器。MTI滤波器的设计目标就是设计一组合适的滤波器系数，使其有效地抑制杂波，并保证目标信号能无损伤地通过。MTI滤波器地优化设计方法主要有特征矢量法和零点分配法。

特征矢量法

特征矢量法是以平均改善因子最大为准则的杂波抑制方法。

通常假设杂波具有高斯型功率谱，谱中心为$f_0$，谱宽为${\delta }_f$，谱密度函数为$C(f)=\frac{1}{2\pi {\delta }_f}{e}^{-\frac{(f-f_0{)}^2}{2{\delta }_f^2}}$

根据维纳滤波理论，如果杂波是平稳随机过程，其功率谱与自相关函数是傅里叶变换对的关系。所以，杂波自相关函数$r_c(m,n)$为其功率谱$C(f)$的傅里叶逆变换

$r_c(m,n)={\int }_{-\infty }^{\infty }C(f)e^{j2\pi f(t_m-t_n)}={\int }_{-\infty }^{\infty }\frac{1}{2\pi {\delta }_f}{e}^{-\frac{(f-f_0{)}^2}{2{\delta }_f^2}}{e}^{j2\pi f(t_m-t_n)}df$

利用积分公式

$\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{-a{x}^2}{e}^{j\xi \tau }dx=\frac{1}{\sqrt{2a}}{e}^{\frac{{\xi }^2}{4a}}$

经推导得到

$r_c(m,n)=e^{-2{\pi }^2{\delta }_f^2{\tau }_{mn}^2}{e}^{j2\pi f_0{\tau }_{mn}}$

式中，${\tau }_{mn}=t_m-t_n$为相关时间。如果杂波谱的中心频率为零，这时

$r_c(m,n)=e^{-2{\pi }^2{\delta }_f^2{\tau }_{mn}^2}$

得到N个脉冲的杂波的自相关矩阵$R_c$为

$$R_c=\left(\begin{array}{cccc}{r}_c(0,0)& r_c(0,1)& \cdots & r_c(0,N-1)\\ r_c(1,0)& r_c(1,1)& \cdots & r_c(1,N-1)\\ ⋮& ⋮& ⋮& \ddots \\ r_c(N-1,0)& r_c(N-1,1)& \cdots & r_c(N-1,N-1)\end{array}\right)$$

对于目标回波信号来说，其多普勒频率是未知的，假设其在区间$[-\frac{B}{2},\frac{B}{2}]$上为均匀分布，且$B>>f_r$，则目标回波信号的多普勒频谱$S(f)$可表示为

$$S(f)=\{\begin{array}{ll}1,& -\frac{B}{2}\le f\le \frac{B}{2}\\ 0,& 其他\end{array}$$

目标的自相关函数为

$$r_s(m,n)=\frac{1}{B}{\int }_{-B/2}^{B/2}{e}^{j2\pi f{\tau }_{mn}}df=\frac{1}{j2\pi B{\tau }_{mn}}[e^{j2\pi B{\tau }_{mn}/2}-e^{-j2\pi B{\tau }_{mn}/2}]=\frac{sin(\pi B{\tau }_{mn})}{\pi B{\tau }_{mn}}=\{\begin{array}{ll}1,& m=n\\ 0,& m\ne n\end{array}$$

假设N脉冲MTI输入端的杂波数据和目标数据分别为

$C=[c(t_1)，c(t_2)，\dots \dots ，c(t_N){]}^T$

$S=[s(t_1)，s(t_2)，\dots \dots ，s(t_N){]}^T$

那么MTI输出端的杂波功率和信号功率分别为

$C_o=E[|w^{H}C{|}^2]=C_i{w}^H{R}_{c}w$

$S_o=E[|w^{H}C{|}^2]=S_i{w}^H{R}_{s}w$

式中，$C_i$和$S_i$分别表示MTI滤波器输入端的杂波功率和信号功率，$w$为FIR滤波器权系数矢量。根据MTI滤波器的改善因子的定义

$I=\frac{S_o/C_o}{S_i/C_i}=\frac{S_o}{S_i}\frac{C_i}{C_o}=\frac{S_i{w}^H{R}_{s}w}{S_i}\frac{C_i}{C_i{w}^H{R}_{c}w}=\frac{w^H{R}_{s}w}{w^H{R}_{c}w}$

由$r_s(m,n)$知，$R_s$为单位阵，因此

$I=\frac{w^{H}w}{w^H{R}_{c}w}$

$R_c$的特征方程为

$R_c{w}_n={\lambda }_n{w}_n，n=0,1,\dots \dots ,N$

式中，$w_n$为特征值${\lambda }_n$所对应的特征向量。其中

${\lambda }_0\le {\lambda }_1\le \dots \dots \le {\lambda }_n$

在$R_c$的特征值中，大特征值所对应的特征向量张成的子空间为信号子空间，杂波的主要分量位于这个子空间；小特征值所对应的特征向量张成的子空间为噪声子空间。因为噪声子空间与信号子空间是正交的，所以最小特征值${\lambda }_0$所对应的特征向量$w_0$被取为MTI滤波器的权系数向量，这就可以在最大程度上抑制杂波分量，使改善椅子最大。

这种利用杂波自相关矩阵的特征分解，用其最小特征值所对应的特征向量设计MTI滤波器的方法称为特征矢量法。这样设计的滤波器可以得到良好的杂波抑制性能，应用较广泛。

如果存在两种或两种以上的杂波，如地杂波和云雨杂波，两种杂波的谱中心可能分别位于不同频率轴上不同位置，对于多个高斯谱的混合杂波，其功率谱是它们各自功率谱之和，其自相关函数也由对应的多杂波分量之和构成。可以用特征矢量法设计具有两个凹口的滤波器，同时在两种杂波谱中心形成两个不同的凹口。

代码：

firjie = 4; %滤波器阶数
fr = 100; %平均脉冲重复频率
bianT = [27,28,29]; %变T码
fd = 0; %杂波谱中心
df = 0.64; %杂波谱宽

bianT = bianT/fr/mean(bianT);
bianT_num = length(bianT);
bianT = repmat(bianT,1,firjie);
T = zeros(1,firjie);
f = [-1*fr:0.1:4*fr]; %频率范围
hd = zeros(bianT_num,length(f));
Rc = zeros(firjie,firjie); %杂波自相关矩阵
for k=1:bianT_num
    for m = 2:firjie
        T(m) = sum(bianT(k:k+m-2));
    end
    for m = 1:firjie
        Rc(m,:) = exp(((fd+1i*2*pi*df^2*(T(m)-T)).^2-fd^2)/(2*df.^2));
    end
[V,D] = eig(Rc);
w = V(:,1);
ww(k,:) = w;
end
hd = ww* exp(-1i*2*pi*T'*f); %计算频率响应
Hd = 20*log10(abs(hd));
figure;
plot(f./fr,Hd);

运行结果：

零点分配法

零点分配法是在设计带阻滤波器时，在凹口处设置频率响应零点的一种方法。在自适应杂波抑制的应用环境下，需要的理想滤波器是在杂波分布的频点处，使杂波得到最大限度地抑制，而在其他频率点具有最大平坦幅度。

对于N阶FIR滤波器，滤波器的权系数为${\omega }_i(i=0,1,2,\dots \dots ，N)$，滤波器的频率响应函数为

$H(f)={\sum }_{i=0}^N{\omega }_i{e}^{-j2\pi f{T}_i}$

式中，$T_0=0,T_i={\sum }_{k=1}^i{T}_{r,k}(i=1,2,\dots \dots ,N)$，$T_{r,k}$为每个脉冲之间的时间间隔。

将$H(f)$在$f=f_0$处展开成泰勒级数

$H(f)=H(f_0)+H^{\prime }(f_0)(f-f_0)+\frac{H^{\prime \prime }(f_0)}{2!}(f-f_0{)}^2+\dots \dots$

式中$H^{(k)}(f_0)=(-j2\pi {)}^k{\sum }_{i=0}^N{T}_i^k{\omega }_i{e}^{-j2\pi f_0{T}_i}$

要在$f=f_0$处设计带阻滤波器，则必须使泰勒级数展开式中$(f-f_0{)}^k$的系数为0，即$H^{(k)}(f_0)=0$。这样就产生了N个关于${\omega }_i(i=0,1,2,\dots \dots ,N)$的齐次线性方程

式中$H^{(k)}(f_0)=(-j2\pi {)}^k{\sum }_{i=0}^N{T}_i^k{\omega }_i{e}^{-j2\pi f_0{T}_i}$

要在$f=f_0$处设计带阻滤波器，则必须使泰勒级数展开式$(f-f_0{)}^k$的系数为0，即$H^{(k)}(f_0)=0$。这样就产生了N个关于${\omega }_i(i=0,1,2,\dots \dots ,N)$的齐次线性方程

${\sum }_{i=0}^N{T}_i^k{\omega }_i{e}^{-j2\pi f_0{T}_i}=0，k=0,1,\dots \dots ,N-1$

其中，$T_i^0=1$，${\omega }_0$为一个常数，通常设置为1。将上式写成矩阵形式

$A\omega =-{\omega }_{0}U$

其中，$w=[w_1，w_2，\dots \dots ，w_N{]}^T$，$U=[1，0，\dots \dots ，0{]}^T$，

$$A=\left(\begin{array}{cccc}{e}^{-j2\pi f_0{T}_1}& e^{-j2\pi f_0{T}_2}& \cdots & e^{-j2\pi f_0{T}_N}\\ T_1{e}^{-j2\pi f_0{T}_1}& T_2{e}^{-j2\pi f_0{T}_2}& \cdots & T_N{e}^{-j2\pi f_0{T}_N}\\ ⋮& ⋮& ⋮& \ddots \\ T_1^{N-1}{e}^{-j2\pi f_0{T}_1}& T_2^{N-1}{e}^{-j2\pi f_0{T}_2}& \cdots & T_N^{N-1}{e}^{-j2\pi f_0{T}_N}\end{array}\right)$$

可以使用Gauss-Jordan法求解这一方程，得到当前时刻的滤波器权系数。

当阻带处于零频时，所设计的滤波器为零频处最大平坦阻带滤波器，此时A为Vandermonde矩阵

$$A=\left(\begin{array}{cccc}1& 1& \cdots & 1\\ T_1& T_2& \cdots & T_N\\ ⋮& ⋮& ⋮& \ddots \\ T_1^{N-1}& T_2^{N-1}& \cdots & T_N^{N-1}\end{array}\right)$$

在单零点时，滤波器阻带凹口较窄。为此，可以在阻带内多设置几个零点，使其阻带拓宽。当然，滤波器的长度应大于零点的个数。

代码：

mtijie = 4; %MTI滤波器的长度
fr = 100; %平均脉冲重复频率
bianT = [27,28,29]; %变T码
fz = 0; %零点位置
amtijie = mtijie;
bianT = bianT/fr/mean(bianT);
bianT_num = length(bianT);
f = (-fr:0.1:4*fr); %频率范围
coenum = amtijie - 1; %需计算的权值个数=FIR长度-1；第一个权值为1，不需要计算
tao_T = repmat(bianT,1,coenum);
w = zeros(1,amtijie); %权系数
for m = 1:bianT_num
    w(1) = 1; %第一个权系数为常数
    U = [1;zeros(coenum-1,1)]; %U矩阵
    tao = tao_T(m:m+coenum);
    Ti(1) = 0;
    for n = 1:coenum
        Ti(n+1) = sum(tao(1:n));
    end
    for n = 1:coenum
        A(n,1:coenum)=(Ti(2:end).^(n-1)).*exp(-1i*2*pi*fz.*Ti(2:end));
    end
    w(2:end) = -w(1)*inv(A)*U;
    w = w./max(abs(w)); %归一化
    ww(m,:) = w;
end
hd = ww * exp(-1i*2*pi*Ti'*f); %计算频率响应
Hd = 20*log10(abs(hd));
figure;
plot(f./fr,Hd);
ylim([-100,10]);

运行结果：

参考文献：

1. 《现代雷达系统分析与设计》，陈伯孝