信号处理工具箱

一、信号的表示 在信号处理中大多数信号是需要事先输入时间向量的，对于只有一个输入和一个输出的信号，MATLAB 将通过向量的形式来表示它。假设输入为时间 t，输出信号为 y。取 t=t1:p:t2，其中 t1 表示时间信号的起始时间，t2 表示时间信号的终止时间，p 为时间间隔，那么输出信号 y=sin(t)可以由时间向量 t 和 t 向量在 f(t)对应时间点上的采样值表示。e.g.y=sin(t)可以由时间向量 t 和向量 y 表示 t=-10:1:10%输入时间向量，按下回车键显示时间向量如下 t=Columns 1 through 11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 Columns 12 through 21 1 2 3 4 5 6 7 8 9 10 y=sin(t)%输入连续正弦信号，按下回车键显示 y 向量 y=Columns 1 through 6 0.5440 -0.4121 -0.9894 -0.6570 0.2794 0.9589 Columns 7 through 12 0.7568 -0.1411 -0.9093 -0.8415 0 0.8415 Columns 13 through 18 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.6570 Columns 19 through 21 0.9894 0.4121 -0.5440 使用绘图命令 plot(t,y)，可以看到由向量 t 和 y 表示的正弦信号如图 1.1 所示。-10-8-6-4-20246810-1-0.8-0.6-0.4-0.200.20.40.60.81-10-8-6-4-20246810-1-0.8-0.6-0.4-0.200.20.40.60.81 可见，其图像有些失真，这是因为 plot 命令是将图中间隔两点用直线连接的，若减小时间间隔 p，将有效的恢复正弦函数原貌。t=-10:0.01:10;y=sin(t);plot(t,y)二、信号的生成 1、正弦波和余弦波 在 MATLAB 中利用函数 sin 和 cos 可以生成所需要的正弦波或余弦波。e.g.生成一个信号持续时长 10s，频率为 250Hz，幅度为 0.75，初始相位为 40的余弦波，并画出其波形图。（1）问题分析：根据采样定理（采样速率大于等于模拟信号的最高频率的 2 倍，模拟信号可以由采样序列构成的时间离散信号无失真的表达），采样率至少为信号最高频率的两倍。采样频率必须大于 500 次/秒，为了产生光滑的曲线，本例中取采样频率为 3000 次/秒，信号持续时间为 10 秒，那么采样点数为 103000=30000。（2）MATLAB 命令生成。Fs=3000;%采样频率 t=1/Fs:1/Fs:10;%信号的持续时间 f=250;%余弦波频率 A=0.75;%信号幅度 Ip=40/180*pi;%初始相位 y=A*cos(2*pi*f.*t+Ip);%余弦波计算 plot(t(1:100),y(1:100)%画出余弦波前 100 个采样值（3）由余弦波前 100 点采样值绘出的图形 00.0050.010.0150.020.0250.030.035-0.8-0.6-0.4-0.200.20.40.60.8 2、周期方波和锯齿波 square 命令生成方波，sawtooth 命令生成三角波，也称锯齿波。它们的调用格式如下：square(T)：产生一周期为 2，幅值为 1 的方波，采样频率由向量 T 指定；square(T,DUTY)：产生一个给定占空比，周期为 2，幅值为 1 的方波，占空比是 1100 之间的数，如果占空比是 30，表示一个方波的周期内正电平占 30%；sawtooth(T)：产生周期为 2，幅值为 1 的三角波，采样时刻由向量 T 指定；sawtooth(T,WIDTH)：产生三角波，WIDTH 指定最大值出现的地方，其取值在 0 到 1 之间。当 T 由 0 增大到 WIDTH*2 时，函数值由-1 增大到 1，当 T 由 WIDTH*2 增大到 2 时，函数值由 1 减小到-1。3、周期 sinc 函数 周期 sinc 函数在 MATLAB 中用 diric 命令实现，其又称为 Dirichlet 函数。Dirichlet 函数定义：d(x)=sin(N*x/2)./(N*sin(x/2)。diric()函数的调用格式：Y=diric(X,N)，其返回的是一个大小与 X 相同的矩阵，其元素为Dirichlet 函数值。N 必须为正整数，改函数将 0 到 2 等间隔的分成 N 等份。e.g.生成 sinc 函数。（1）MATLAB 程序实现。x=0:0.03:3*pi;y1=diric(x,5);y2=diric(x,9);subplot(121);plot(x,y1);xlabel(x);title(N=5);subplot(122)plot(x,y2);xlabel(x);title(N=9);（2）程序运行结果 0510-0.4-0.200.20.40.60.81xN=50510-0.4-0.200.20.40.60.81xN=9 sinc 函数不同 N 值下的波形图 4、高斯调整正弦脉冲 Gauspuls 是 MATLAB 信号处理工具箱提供的信号发生函数，其调用格式如下：YI=gauspuls(T,FC,BW)：函数返回最大幅值为 1 的高斯函数调幅的正弦波的采样，其中心频率为FC，相对带宽为BW，时间由数组T给定。BW的值必须大于0。默认情况下，FC=1000Hz，BW=0.5。YI=gauspuls(T,FC,BW,BWR)：BWR 指定可选的频带边缘处的参考水平，以相对于正常信号峰值下降了-BWR（单位为 dB）为边界的频带，其相对带宽为 100*BW%。默认情况下，BWR的值为-6dB。其他参数设置同上。BWR 的值为负值。TC=gauspuls(cutoff,FC,BW,BWR,TPE)：返回包络相对包络峰值下降 TPE（单位为 dB）时的时间 TC。默认情况下，TPE 的值是-60dB。其他参数设置同上。TPE 的值必须是负值。e.g.生成一个中心频率为 200kHz 的高斯调幅正弦脉冲，其相对带宽为 0.6，采样率为20MHz，信号在包络相对于峰值下降 30dB 时截断。（1）MATLAB 程序实现 tc=gauspuls(cutoff,200e3,0.6,-30);t=-tc:1e-7:tc;yi,yq,ye=gauspuls(t,200e3,0.6);plot(t,yi,t,yq,t,ye)（2）程序运行结果-1-0.8-0.6-0.4-0.200.20.40.60.81x 10-5-1-0.8-0.6-0.4-0.200.20.40.60.81 5、扫频信号 利用 MATLAB 的信号处理工具箱中的 chirp 函数可以获得在设定频率范围内的按照设定方式进行的扫频信号，调用格式如下。Y=chirp(T,R0,T1,F1)：产生一个频率随时间线性变化信号的采样，其时间轴的设置由数组 T定义，时刻 0 的瞬时频率为 F0；时刻 T1 的瞬时频率为 F1。默认情况下，F0=0Hz，T1=1，F1=100Hz。Y=chirp(T,R0,T1,F1,method)：method 指定改变扫频的方法。可用的方法有linear（线性调频）、quadratic（二次调频）、logarithmic（对数调频），默认时为linear。Y=chirp(T,R0,T1,F1,method,PHI)：PHI 指定信号的初始相位，默认时 PHI 的值为 0。e.g.以 1000Hz 的采样频率，在 3s 采样时间内，生成一个起始时刻瞬时频率是 10Hz，3s 时瞬时频率为 100Hz 的线性调频信号，并画出其曲线图及光谱图。MATLAB 程序实现如下：fs=1000;t=0:1/fs:3;y=chirp(t,0,1,100);plot(t(1:500),y(1:500);%画出线性调频信号的前六分之一段 t1=0:1/fs:3;%信号光谱图的实现 y=chirp(t1,0,1,100);specgram(y,256,1e3,256,250);%线性调频信号的光谱图 00.050.10.150.20.250.30.350.40.450.5-1-0.8-0.6-0.4-0.200.20.40.60.81TimeFrequency0.511.522.5050100150200250300350400450500 6、单位冲激信号 单位冲激信号是信号系统的基本信号，表示符号(t)，数学定义是：-1)(0,0)(dtttt 单位冲激信号除了原点之外，其他各处都为零，并且信号的总面积为 1，实际就用一个矩形脉冲来代替单位冲激信号。单位冲激信号在 MATLAB 中的实现代码如下：dt=0.01;t=-3:dt:3;%画出-3 3上的波形图 n=length(t);%计算采样点数 n x=zeros(1,n);%生成一维数组，其 n 个元素都为 0 x(1,3/dt+1)=1/dt;%设置原点处的采样值 stairs(t,x);%画出信号图 axis(-3,3,0,150)%设置显示窗口横纵坐标 -3-2-10123050100150 将坐标设置改为 axis(-0.03,0.03,0,150)后可以看到单位冲激信号为窄矩形脉冲近似表示。-0.0 3-0.0 2-0.0 100.0 10.0 20.0 3050100150 7、单位序列 单位序列是一个典型的离散信号，数学表达0010)(kkk 。只有当 k=0 时，函数值才为 1，其余全为零，在 MATLAB 中实现较为简单，下面将给出平移信号(k+k0)的 MATLAB实现代码，改信号是指在 k=-k0时函数值为 1，其余时刻函数值为 0。e.g 画出单位序列 (k+k0)在k1 k2上的波形图，假设 k1=-5，k2=6，k0=3。MATLAB 实现代码如下：k1=-5;k2=6;k0=3;k=k1:k2;n=length(k);x=zeros(1,n);x(1,-k0-k1+1)=1;stem(k,x,filled);%绘制波形 axis(k1,k2,-0.3,1.3);xlabel(k);title(单位序列)-5-4-3-2-10123456-0.200.20.40.60.811.2k单 位 序 列 8、均匀分布的随机序列 rand 函数可以生成在0 1区间上均匀分布的随机数序列，rand 函数的一般调用格式为：Y=rand(M,N)，其生成 M 行 N 列的随机数矩阵。e.g 产生一个在-3 5上均匀分布的 6 行 5 列的随机数矩阵，并画出其随机数发生频率分布图。MATLAB 中实现命令如下：M=6;N=5;a=-3;b=5;A=rand(M,N)*(b-a)+a%生成-3 5上均匀分布的随机数矩阵，共 5 行 6 列 A=%生成的随机数矩阵 3.5178 -0.7720 4.6573 3.3377 2.4299 4.2463 1.3751 0.8830 4.6759 3.0619 -1.9841 4.6601 3.4022 2.2459 2.9451 4.3070 4.7191 -1.8649 -2.7143 0.1378 2.0589 -1.7391 0.3741 3.7930 2.2438 -2.2197 4.7647 4.3259 4.4719 -1.6305 A=rand(1,100000)*(b-a)+a;%检验生成随机数的范围 min(A)%最小值 ans=-2.9999 max(A)%最大值 ans=4.9998 x=-4:0.1:6;%画出 A 中随机数的发生频率分布 hist(A,x)-6-4-2024680200400600800100012001400 9、高斯分布随机序列 在 MATLAB 中，除了均匀分布的随机序列外，常用的还有标准正态分布的随机序列，该序列可以用 randn 函数生成，randn 函数的调用格式：Y=randn(M,N)，将生成一个 M 行 N 列的均值方差为 1 的标准正态分布的随机数序列。e.g.产生一个 10000 个均值为 100、方差为 4.6 的正态分布的随机数序列，并画出其随机数发生频率分布图。MATLAB 中实现命令如下：M=100;%设置均值 D=4.6;%设置方差 Y=M+sqrt(D)*randn(1,10000);%生成随机数序列 M1=mean(Y)%检验生成序列均值 M1=99.9955 D1=var(Y)%检验生成序列方差 D1=4.6077 x=80:0.1:120;hist(Y,x)%绘制 Y 中随机数发生频率分布图 7580859095100105110115120125050100150200250 三、随机信号处理和谱分析 1、随机信号互相关函数估计 xcorr 函数是随机信号互相关估计函数，其调用格式如下。c=xcorr(x,y,maxlags,option)：上式中，x，y 为两个独立的随机信号序列，长度都为 N 的向量（N1），如果 x，y 长度不一致，则短的一个用 0 补齐，使得两个信号长度一样；c 为 x，y 的互相关函数估计序列；maxlags 为 x 与 y 之间的最大延迟，函数返回值 c 的长度为2maxlags+1。默认状态下，函数返回值 c 的长度为 2N-1。option指定互相关的归一化选项，它可以是：biased：计算互相关函数的有偏互相关估计；unbiased：计算互相关函数的无偏互相关估计；coeff：系列归一化，使零延迟的自相关为 1；none：默认状态，函数执行非归一化计算相关。e.g.已知两个信号 x(t)=sin(2 ft)，y(t)=kcos(2 ft+)式中，f 为 10Hz，k 为 2，为/6。求这两个信号的自相关函数 Rx()、Ry()以及互相关函数 Rxy()。（1）MATLAB 程序实现如下：Fs=1000;N=1000;n=0:N-1;t=n/Fs;Lag=200;f=10;k=2;w=pi/6;x=sin(2*pi*f*t);%信号 x(t)y=k*cos(2*pi*f*t+w);%信号 y(t)cx,lagsx=xcorr(x,Lag,unbiased);%求 x 的自相关函数 cy,lagsy=xcorr(x,Lag,unbiased);%求 y 的自相关函数 c,lags=xcorr(x,y,Lag,unbiased);%求 xy 互相关函数 subplot(311);%画出信号 x(t)的自相关函数 plot(lagsx/Fs,cx,r);xlabel(t);ylabel(Rx(t);title(信号 x 自相关函数);subplot(312)%画出信号 y(t)的自相关函数 plot(lagsy/Fs,cy,b);xlabel(t);ylabel(Ry(t);title(信号 y 自相关函数);subplot(313);%画出互相关函数 plot(lags/Fs,c,r);xlabel(t);ylabel(Rxy(t);title(互相关函数);grid;（2）程序运行结果 -0.2-0.15-0.1-0.0500.050.10.150.2-101tRx(t)信 号 x自 相 关 函 数-0.2-0.15-0.1-0.0500.050.10.150.2-101tRy(t)信 号 y自 相 关 函 数-0.2-0.15-0.1-0.0500.050.10.150.2-202tRxy(t)互 相 关 函 数 2、互协方差函数估计 xcov 函数是互协方差估计函数，其调用格式如下。c,lags=xcov(x,y,maxlags,option)：上式中，x，y 为两个独立的随机信号序列，长度都为 N的向量（N1）；v 为 x，y 的互协方差序列；maxlags 为 x 与 y 之间的最大延迟，函数返回值 c 的长度为 2N-1。option指定互协方差的归一化选项，它可以是：biased：计算互协方差函数的有偏互相关估计；unbiased：计算互协方差函数的无偏互相关估计；coeff：系列归一化，使零延迟的自相关为 1；none：默认状态，函数执行非归一化计算相关。e.g.估计一个正态分布白噪声信号 x 的自协方差 cx(n)，假设最大延迟为 30。（1）MATLAB 程序实现如下：x=randn(1,1000);%生成白噪声 cov_x,lags=xcov(x,30,coeff);%计算协方差 stem(lags,cov_x)%图像输出（2）程序运行结果-30-20-100102030-0.200.20.40.60.811.2 3、谱分析函数 psd 对于随机信号来说，虽然它没有确定的解析表达式，但其相关函数却是确定的，如果信号是平稳的，那么对相关函数的傅里叶变换就是它的功率谱密度函数，及功率谱，它反映了单位频带内随机信号功率的大小。MATLAB 中提供的谱分析函数最常用的是函数 psd 和 pwelch求功率谱。功率谱密度估计函数（psd 函数）的调用格式如下。Pxx=psd(X,NFFT,Fs,WINDOW)：返回的是信号向量 X 的功率谱密度估计，用的是 Welch 平均周期图法。X 被分成重叠的几段，每段都经过了去趋势，再经过由参数 WINDOW 指明的窗函数加窗，然后补零至 NFFT 点长。各段 NFFT 点 DFT 的幅值的平方的平均值即为 Pxx。Pxx 的长度：当 NFFT 为偶数时，其值为 NFFT/2+1；当 NFFT 为奇数时，其值为（NFFT+1）/2；当 NFFT 为复数时，为 NFFT。当 WINDOW 为一数值 n 时，则采用 n 点长的 Hanning窗加窗。Pxx,F=psd(X,NFFT,Fs,WINDOW,NOVERLAP)：返回一由频率点组成的向量 F，Pxx 为这些点上的估值，X 在分段时相邻两段有 NOVERLAP 点重叠。其他参数同上。Pxx,Pxxc,F=psd(X,NFFT,Fs,WINDOW,NOVERLAP,P)：返回 Pxx 的 P*100%置信区间 Pxxc，其中 P 在 0 到 1 间取值。其他参数同上。Psd(X,NFFT,Fs,WINDOW,NOVERLAP,P,DFLAG)：参数 DFLAG 可选：linearmeannone。DFLAG 指明了在对各段加窗后进行趋势去除时使用的方式。默认或为空矩阵的情况下，NFFT 为 256，若 x 的长度小于 256 时为 X 的长度；NOVERLAP 为 0；WINDOW为 Hanning，数值表述为 NFFT；Fs 为 2；P 为 0.95；DFLAG 为none。其他参数同上。e.g.采用采样频率为 1000Hz、长度为 1024 点、相邻两段重叠点数为 128，窗函数为默认值的 Welch 方法对信号 x(t)=sin2 f1t+ksin2 f2t+n(t)进行功率谱估计，其中 n(t)正态分布白噪声，f1=100Hz，f2=200Hz，k=2。（1）MATLAB 程序实现如下：fs=1000;%设置采样频率 t=0:1/fs:1;f1=100;f2=200;k=2;x=sin(2*pi*f1*t)+k*sin(2*pi*f2*t)+randn(1,length(t);%信号 x(t)p,f=psd(x,1024,1000,128);%功率谱密度估计 plot(f,10*log10(p/(1024/2);%绘制功率谱密度图 xlabel(freq Hz);ylabel(PSD);title(功率谱估计)（2）程序运行结果 050100150200250300350400450500-60-50-40-30-20-10010freq HzPSD功 率 谱 估 计 4、谱分析函数 pwelch pwelch 函数的调用格式如下。Pxx,w=pwelch(x)：该函数用 Welch 方法估计一个输入信号向量 x 的功率谱密度 Pxx。其向量 x 被分割成 8 段，每一段有 50%的重叠，函数将忽略没有包含在 8 段中的剩余的 x 中的数据，并且这分割后的每一段都用汉明窗进行加窗，窗函数的长度和每一段的长度一样。当x 为实数时，产生单边的 PSD，当 x 是复数时，产生双边的 PSD。一般来说，FFT 的长度和输入 x 的值决定了 Pxx 的长度和归一化频率 w 的范围，系统默认 FFT 的长度 N 为 256 和 2的整数次幂中大于分段长度的最近的数。具体规定为，当输入x是实数时，Pxx的长度为（N/2）+1，对应的归一化频率的范围为0,；当输入 x 是复数时，Pxx 的长度为 N，对应的归一化频率范围为0,2)。Pxx,w=pwelch(x,window)：如果设定 window 是一个正整数，那么这个数代表 Hamming 窗的长度；如果设定 window 为一个向量，那么这个向量代表窗函数的权系数。在这种调用格式中，输入向量 x 被分割成每段重叠 50%的整数段，每段的长度和窗函数的长度相同，没有包含在任何一段中的剩余的 x 中的数据将被忽略。如果指定 window 为一个向量，则信号数据被分割成 8 段，并在每一段上加 Hamming 窗。Pxx,w=pwelch(x,window,noverlap)：改调用格式指定 x 分割后每一段的长度为 window，noverlap 指定每段重叠的信号点数，noverlap 必须小于被确定的窗口长度，在默认情况下，x被分割后的每段有 50%重叠。Pxx,w=pwelch(x,window,noverlap,nfft)：整数 nfft 指定 FFT 的长度，如果 nfft 指定为一个空向量，则 nfft 取前面调用格式中的 N。nfft 和 x 决定了 Pxx 的长度和 w 的频率范围，具体规定为：当输入 x 为实数、nfft 为偶数时，Pxx 的长度为（nfft/2+1），w 的范围为0,；当输入 x 为实数、nfft 为奇数时，Pxx 的长度为（nfft+1）/2，w 的范围为0,；当输入 x 为复数、nfft 为偶数或奇数时，Pxx 的长度为 nfft，w 的范围为0,2。Pxx,f=pwelch(x,window,noverlap,nfft,fs)：整数 fs 为采样频率，如果定义 fs 为空间向量，则采样频率默认为 1Hz。nfft 和 x 决定 Pxx 的长度和 f 的频率范围，具体规定为：当输入 x 为实数、nfft 为偶数时，Pxx 的长度为（nfft/2+1），f 的范围为0,fs/2；当输入 x 为实数、nfft为奇数时，Pxx 的长度为（nfft+1）/2，f 的范围为0,fs/2)；当输入 x 为复数、nfft 为偶数或奇数时，Pxx 的长度为 nfft，f 的范围为0,fs)。.=pwelch(x,window,noverlap,.,range)：当 x 是一个实数的时候，这种调用格式非常有用，它确定 f 或 w 的频率取值范围。字符串 range 可以取 twosided 和 onesided。twosided计算双边 PSD；onesided计算单边 PSD。pwelch(.)：该命令在当前 Figure 窗口中绘制出功率谱密度曲线，其单位为 dB/Hz。e.g.估计信号 x(t)=sin2 f1t+ksin2 f2t+n(t)的功率谱密度，显示双边 PSD。采样频率为1000Hz、窗口长度为 60 点、相邻两段重叠点数为 40，FFT 长度为默认值，其中 n(t)正态分布白噪声，f1=100Hz，f2=200Hz，k=2。（1）MATLAB 程序实现如下：fs=1000;%设置采样频率 t=0:1/fs:1;f1=100;f2=200;k=2;x=sin(2*pi*f1*t)+k*sin(2*pi*f2*t)+randn(1,length(t);%信号 x(t)pwelch(x,60,40,fs,twosided)（2）运行程序结果 0100200300400500600700800900-32-30-28-26-24-22-20-18-16-14-12Frequency(Hz)Power/frequency(dB/Hz)Welch Power Spectral Density Estimate 四、模拟滤波器设计 对特定频率的频点或该频点以外的频率进行有效滤除的系统，就是滤波器，换句话说，滤波器可以被用来区分不同频率的信号，实现各种模拟信号的处理过程。介绍内容：模拟滤波器的设计参数、常用的两种模拟滤波器（巴特沃斯滤波器和切比雪夫滤波器）、模拟滤波器的相互转换。1、滤波器的设计参数 滤波器有四个常用的参数，它们分别是：Fp：通带截止频率，单位 Hz；Fs：阻带起始频率，单位 Hz；Rs：阻带内最小衰减，又称阻带波纹，单位 dB；Rp：通带内允许的最大衰减，又称通带波纹，单位 dB。由此，假设采样频率为 fN，那么将得到另外两个归一化角频率参数：p:通带截止角频率，单位 rad/s，计算公式为p=fp/(fN/2)；s:阻带起始角频率，单位 rad/s，计算公式为s=fs/(fN/2)。通过这些参数就可以衡量一个滤波器的性能了，当给定这些参数时，滤波器的设计参数就给定了，用户依照这些参数进行设计。2、巴特沃斯滤波器 巴特沃斯模拟低通滤波器原型的平方幅频响应函数为NcAjH222)(11)()(，c是低通滤波器的截止频率，N 为滤波器的阶数。N 越大，通带和阻带的近似性越好，过渡带也越陡，特性越接近于矩形。巴特沃思滤波器有以下特点：=c时，A(c2)/A(0)=1/2，幅度衰减21/，相当于 3dB 衰减点；在/c1，即在过渡带和阻带中，A（2）单调减小，因为/c1，所以 A（2）快速下降。在MATLAB中butterworth模拟滤波器函数调用格式为Z,P,K=buttap(N)：N 为butterworth滤波器的阶数，Z，P，K 分别为滤波器的零点、极点、增益。该函数返回 N 阶低通模拟滤波器原型的极点和增益。e.g.画出 5 阶和 10 阶 butterworth低通滤波器的平方幅频响应曲线。（1）MATLAB程序实现如下：n=0:0.01:2;N1=5;%设置阶数 N2=10;z1,p1,k1=buttap(N1);%模拟低通滤波器函数 z2,p2,k2=buttap(N2);b1,a1=zp2tf(z1,p1,k1);%零极点模型转化为传递函数模型 b2,a2=zp2tf(z2,p2,k2);H1,w1=freqs(b1,a1,n);H2,w2=freqs(b2,a2,n);magH1=(abs(H1).2;magH2=(abs(H2).2;plot(w1,magH1,w2,magH2);axis(0 2 0 1);xlabel(w/wc);ylabel(|H(jw)|2)（2）程序运行结果 00.20.40.60.811.21.41.61.8200.10.20.30.40.50.60.70.80.91w/wc|H(jw)|2 3、切比雪夫型滤波器 巴特沃斯滤波器的频率特性无论在通带与阻带都随频率而单调变化，因而如果在通带边缘满足指标，则在通带内肯定满足，也就是会超过指标的要求，因而并不经济。所以，更有效的办法是将指标的精度要求均匀地分布在通带内，或均匀分布在阻带内，或同时均匀地分布在通带与阻带内，这样在同样通带、阻带性能要求下，就可设计出阶数较低的滤波器。这种滤波器均匀分布的办法可通过选择具有等波纹特性的逼近函数来完成。切比雪夫滤波器的幅度特性就在一个频带中（通带或阻带）具有这种等波纹特性。切比雪夫型滤波器的平方幅频响应函数为：)(11)()(2222cNCAjH。式中，为小于 1 的正数，表示通带内幅频波纹情况；c是截止频率，N 为多项式)(cNC的阶数，其中1),(coshcosh(1),(coscos()x(11xxNxxNCN。MATLAB 中切比雪夫型滤波器函数为：Z,P,K=cheb1ap(N,Rp)：N 为阶数，Z、P、K 分别为滤波器的零点、极点、增益；Rp 为通带波纹。e.g.画出 5 阶切比雪夫型模拟低通滤波器原型的平方幅频响应曲线。（1）MATLAB 程序实现如下：n=0:0.01:2;N=5;Rp=0.5;%设置通带波纹 z,p,k=cheb1ap(N,Rp);%切比雪夫型模拟滤波器原型 b,a=zp2tf(z,p,k);H,w=freqs(b,a,n);magH=(abs(H).2;plot(w,magH);axis(0 2 0 1);xlabel(w/wc);ylabel(|H(jw)|2)title(N=5);（2）程序运行结果 00.20.40.60.811.21.41.61.8200.10.20.30.40.50.60.70.80.91w/wc|H(jw)|2N=5 4、切比雪夫型滤波器 切比雪夫型滤波器的平方幅频响应函数为12222)(11)()(cNCAjH。式中，为小于 1 的正数，表示阻带内幅频波纹情况；c是截止频率，N 为多项式)(cNC的阶数，其中1),(coshcosh(1|),(coscos()(11xxNxxNxCN。MATLAB 中切比雪夫型滤波器函数为：Z,P,K=cheb2ap(N,Rs)：N 为阶数，Z、P、K 分别为滤波器的零点、极点、增益；Rs 为阻带波纹。e.g.画出 5 阶切比雪夫型模拟低通滤波器原型的平方幅频响应曲线。（1）MATLAB 程序实现如下：n=0:0.01:2;N=5;Rs=8;%设置通带波纹 z,p,k=cheb2ap(N,Rs);%切比雪夫型模拟滤波器 b,a=zp2tf(z,p,k);%零极点模型转化为传递函数模型 H,w=freqs(b,a,n);%使用基于 FFT 的算法计算模拟滤波器的频率响应 magH=(abs(H).2;plot(w,magH);axis(0 2 0 1.1);xlabel(w/wc);ylabel(|H(jw)|2);title(N=5);（2）程序运行结果 00.20.40.60.811.21.41.61.8200.20.40.60.81w/wc|H(jw)|2N=5 5、模拟滤波器的频域变换 如果想设计其他类型的频率选择性滤波器，如高通、带通、带阻滤波器，可以将一个低通滤波器的频带进行变换，以使其能够表现其他频率选择性滤波器特性。lp2lp 函数，低通滤波器转换为低通滤波器。bt,at=lp2lp(b,a,Wo)：这是函数 lp2lp 的传递函数形式，该函数能改变低通滤波器的截止频率，Wo 为低通滤波器所要实现的截止频率。At,Bt,Ct,Dt=lp2lp(A,B,C,D,Wo)：这是函数 lp2lp 的连续时间状态空间形式。lp2hp 函数：低通滤波器转换为高通滤波器。bt,at=lp2hp(b,a,Wo)：这是该函数的传递函数形式，参数 b、a 是传递函数多项式的分子分母系数，Wo 为转换后高通滤波器的截止频率，函数将模拟低通滤波器转换为指定带宽和中心频率的高通滤波器。At,Bt,Ct,Dt,=lp2hp(A,B,C,D,Wo)：将连续时间状态空间低通滤波器原型转换为一个截止频率为 Wo 的高通滤波器。lp2bp 函数，低通滤波器转换为带通滤波器。bt,at=lp2bp(b,a,Wo,Bw)：这是该函数的传递函数形式，函数将模拟低通滤波器转换为指定带宽和中心频率的带通滤波器，参数 b、a 是传递函数多项式的分子分母系数，Wo 为转换后带通滤波器中心频率，Bw 为转换后带通滤波器带宽。At,Bt,Ct,Dt,=lp2bp(A,B,C,D,Wo,Bw)：将连续时间状态空间低通滤波器原型转换为一个中心频率为 Wo，带宽为 Bw 的带通滤波器。lp2bs 函数，低通滤波器转换为带阻滤波器。bt,at=lp2bs(b,a,Wo,Bw)：这是该函数的传递函数形式，函数将模拟低通滤波器转换为指定带宽和中心频率的带阻滤波器，参数 b、a 是传递函数多项式的分子分母系数，Wo 为转换后带阻滤波器中心频率，Bw 为转换后带阻滤波器带宽。At,Bt,Ct,Dt,=lp2bs(A,B,C,D,Wo,Bw)：将连续时间状态空间低通滤波器原型转换为一个中心频率为 Wo，带宽为 Bw 的带阻滤波器。e.g.设计一个 10 阶 Chebyshev型低通滤波器，Rp=5dB，截止频率 Wo=40Hz，然后将其转换为中心频率 Wo=320Hz、带宽为 Bw=10000Hz 的带通滤波器，并分别画出对应的幅频和相频响应曲线。（1）MATLAB 程序实现如下：n=0:0.01:2;N=10;Rp=5;%设置通带波纹 z,p,k=cheb1ap(N,Rp);%切比雪夫型模拟低通滤波器原型 A1,B1,C1,D1=zp2ss(z,p,k);%零极点模型转化为状态方程模型 Wo1=40*pi;%截止频率 AT1,BT1,CT1,DT1=lp2lp(A1,B1,C1,D1,Wo1);%模拟低通滤波器转换为低通滤波器 num1,den1=ss2tf(AT1,BT1,CT1,DT1);%状态方程模型转化为传递函数模型 figure;freqs(num1,den1);%绘制频率响应曲线 Wo2=320*pi;Bw=10000;%截止频率 AT2,BT2,CT2,DT2=lp2bp(A1,B1,C1,D1,Wo2,Bw);%模拟低通滤波器转换为带通滤波器 num2,den2=ss2tf(AT2,BT2,CT2,DT2);%状态方程模型转化为传递函数模型 figure;freqs(num2,den2);%绘制频率响应曲线 100101102103104-200-1000100200Frequency(rad/s)Phase(degrees)10010110210310410-2010-10100Frequency(rad/s)Magnitude 100101102103104105-200-1000100200Frequency(rad/s)Phase(degrees)10010110210310410510-1510-1010-5100Frequency(rad/s)Magnitude 五、IIR 数字滤波器设计 数字滤波器是数字信号处理的基础，用来对信号进行过滤、检测与参数估计等处理，在通信、图像、语音、雷达等许多领域都有着十分广泛的应用目前数字滤波器的设计有许多现成的高级语言设计程序，但它们都存在设计效率较低，不具有可视图形，不便于修改参数等缺点，而 MATLAB 为数字滤波的研究和应用提供了一个直观、高效、便捷的通道。MATLAB 中自带了许多 IIR 数字滤波器设计函数。1、巴特沃思数字滤波器设计（butter 函数）butter 函数即可以用来设计低通、带通、高通和带阻的模拟巴特沃斯滤波器，也可以用来设计这四类的数字巴特沃思滤波器。数字滤波器设计的 butter 函数调用格式如下：b,a=butter(n,Wn)：函数用来设计一个截止频率为 Wn 的 n 阶低通滤波器。它返回滤波器系数向量 a、b 的长度为 n+1，这些系数按 z 的降幂排列为：nnznazaaznbzbzAzBzH)1(.)2()1()1(.)2()1(b)()()(11 截止频率是滤波器的幅度响应为21/时的频率，归一化截止频率 Wn 取值在0 1之间，这里 1 对应内奎斯特频率。如果 Wn 是一个二元向量，如 Wn=w1 w2；那么函数返回一个带通为 w1ww2、阶数为 2n 的带通数字滤波器。b,a=butter(n,Wn,ftype)：该函数格式设计一个截止频率为 Wn 的高通或者带通数字滤波器，ftype为滤波器类型参数，high为高通滤波器；stop为带阻滤波器。对于带阻滤波器，如果 Wn 是一个二元向量，如 Wn=w1 w2，那么 butter 函数返回一个带通为w1w N=9;Wn=300/500;b,a=butter(N,Wn,high);%巴特沃思高通数字滤波器函数 freqz(b,a,128,1000);%绘制频率响应曲线 title(N=9);%显示阶数（2）程序运行结果 050100150200250300350400450500-1000-5000500Frequency(Hz)Phase(degrees)050100150200250300350400450500-400-300-200-1000Frequency(Hz)Magnitude(dB)N=9 e.g.设计一个 10 阶的巴特沃思带通滤波器，通带频率为 100Hz 到 200Hz，并画出其脉冲响应曲线。（1）MATLAB 程序实现如下 n=5;Wn=100 200/500;b,a=butter(n,Wn);y,t=impz(b,a,101);stem(t,y)（2）程序运行结果 0102030405060708090100-0.25-0.2-0.15-0.1-0.0500.050.10.150.2 2、切比雪夫型数字滤波器设计（cheby1 函数）cheby1 函数即可以用来设计低通、带通、高通和带阻的模拟切比雪夫型滤波器，也可以用来设计这四类的数字切比雪夫型滤波器。它的特性就是通带内等波纹，阻带内单调。cheby1 函数调用格式如下。b,a=cheby1(n,Rp,Wn)：函数用来设计一个截止频率为 Wn，通带波纹为 Rp(dB)的 n 阶切比雪夫低通滤波器。它返回滤波器系数向量 a、b 的长度为 n+1，这些系数按 z 的降幂排列为：nnznazaaznbzbzAzBzH)1(.)2()1()1(.)2()1(b)()()(11 归一化截止频率是滤波器的幅度响应为-Rp(dB)时的频率，对于 cheby1 来说，归一化截止频率 Wn 取值在0 1之间，这里 1 对应内奎斯特频率。如果 Wn 是一个二元向量，如 Wn=w1 w2，那么 cheby1 函数返回一个带通为 w1ww2、阶数为 2n 的带通数字滤波器。b,a=cheby1(n,Rp,Wn,ftype)：该函数格式设计一个截止频率为 Wn、通带波纹为 Rp(dB)的高通或者带通数字滤波器，ftype为滤波器类型参数，high为高通滤波器；stop为带阻滤波器。对于带阻滤波器，如果 Wn 是一个二元向量，如 Wn=w1,w2，那么 cheby1函数返回一个带通为 w1w N=9;Wn=300/500;Rp=0.5;b,a=cheby1(N,Rp,Wn);%切比雪夫型低通数字滤波器函数 freqz(b,a,512,1000);%绘制频率响应曲线 axis(0 500-300 50);title(N=9)%显示阶数（2）程序运行结果 050100150200250300350400450500-1000-5000Frequency(Hz)Phase(degrees)050100150200250300350400450500-300-200-1000Frequency(Hz)Magnitude(dB)N=9 e.g.设计一个 10 阶的切比雪夫型带通滤波器，通带频率为 100Hz 到 200Hz，Rp=0.5dB，并画出其脉冲响应曲线。（1）MATLAB 程序实现如下：n=5;Wn=100 200/500;Rp=0.5;b,a=cheby1(n,Rp,Wn);y,t=impz(b,a,101);stem(t,y)（2）程序运行结果 0102030405060708090100-0.2-0.15-0.1-0.0500.050.10.150.2 3、切比雪夫型数字滤波器设计（cheby2 函数）Cheby2 函数即可以用来设计低通、带通、高通和带阻的模拟切比雪夫型滤波器，也可以用来设计这四类的数字切比雪夫型滤波器。它的特性就是通带内单调，阻带内等波纹。数字滤波器设计的 cheby2 函数调用格式为：b,a=cheby2(n,Rs,Wn)b,a=cheby2(n,Rs,Wn,ftype)z,p,k=cheby2(n,Rs,Wn,s)z,p,k=cheby2(n,Rs,Wn,ftype)A,B,C,D=cheby2(n,Rs,Wn)A,B,C,D=cheby2(n,Rs,Wn,ftype)模拟滤波器设计的 cheby1 函数调用格式为：b,a=cheby2(n,Rs,Wn,s)b,a=cheby2(n,Rs,Wn,ftype,s)z,p,k=cheby2(n,Rs,Wn,s)z,p,k=cheby2(n,Rs,Wn,ftype,s)A,B,C,D=cheby2(n,Rs,Wn,s)A,B,C,D=cheby2(n,Rs,Wn,ftype,s)e.g.设计一个 9 阶的 cheby2 型低通数字滤波器，采样频率为 1000Hz，Rs=20dB，截止频率为 300Hz（归一化频率值为 0.6），并画出滤波器频率响应曲线。（1）MATLAB 程序实现如下：N=9;Wn=300/500;Rs=20;b,a=cheby2(N,Rs,Wn);%切比雪夫型低通数字滤波器函数 freqz(b,a,512,1000);%绘制频率响应曲线 axis(0 500-80 50);title(N=9);%显示阶数（2）程序运行结果 050100150200250300350400450500-400-2000200Frequency(Hz)Phase(degrees)050100150200250300350400450500-50050Frequency(Hz)Magnitude(dB)N=9 e.g.设计一个 10 阶的切比雪夫型带通滤波器，通带频率为 100Hz 到 200Hz，Rs=20dB，并画出其脉冲响应曲线。（1）MATLAB 程序实现如下：n=5;Wn=100 200/500;Rs=20;b,a=cheby2(n,Rs,Wn);y,t=impz(b,a,101);stem(t,y)（2）程序运行结果 0102030405060708090100-0.2-0.15-0.1-0.0500.050.10.150.2 4、椭圆数字滤波器设计（ellip 函数）ellip 函数即可以用来设计低通、带通、高通和带阻的模拟椭圆滤波器，也可以用来设计这四类的数字椭圆滤波器。它在通带内和阻带内都是等波纹的。数字滤波器设计的 ellip 函数调用格式如下。b,a=ellip(n,Rp,Rs,Wn)：函数用来设计一个截止频率为 Wn、通带波纹为 Rp（dB）、阻带波纹为 Rs（dB）的 n 阶椭圆低通滤波器。它返回滤波器系数向量 a、b 的长度为 n+1，这些系数按 z 的降幂排列为：nnznazaaznbzbzAzBzH)1(.)2()1()1(.)2()1(b)()()(11 归一化截止频率是滤波器的幅度响应为-Rp(dB)处的通带边缘频率，对于 ellip 函数来说，归一化截止频率 Wn 取值在0 1之间，这里 1 对应内奎斯特频率。如果 Wn 是一个二元向量，如 Wn=w1 w2，那么 ellip 函数返回一个带通为 w1ww2、阶数为 2n 的带通数字滤波器。b,a=ellip(n,Rp,Rs,Wn,ftype)：该函数格式设计一个截止频率为 Wn、通带波纹为 Rp(dB)、阻带波纹为 Rs(dB)的高通或者带通数字滤波器，ftype为滤波器类型参数，high为高通滤波器；stop为带阻滤波器。对于带阻滤波器，如果 Wn 是一个二元向量，如 Wn=w1 w2，那么 ellip 函数返回一个带通为 w1w N=6;Wn=300/500;Rs=50;Rp=3;b,a=ellip(N,Rp,Rs,Wn);%椭圆数字滤波器函数 freqz(b,a,512,1000);%绘制频率响应曲线 axis(0 500-100 0);title(N=6);%显示阶数（2）程序运行结果 050100150200250300350400450500-600-400-2000200Frequency(Hz)Phase(degrees)050100150200250300350400450500-100-500Frequency(Hz)Magnitude(dB)N=6 E.g.设计一个 20 阶的椭圆带通滤波器，通带频率为 100Hz 到 200Hz，Rp=0.5，Rs=20，并画出其脉冲响应曲线。（1）MATLAB 程序实现如下：N=10;Wn=100 200/500;Rs=20;Rp=0.5;b,a=ellip(N,Rp,Rs,Wn);%椭圆数字滤波器函数 y,t=impz(b,a,101);stem(t,y)（2）程序运行结果 0102030405060708090100-0.2-0.15-0.1-0.0500.050.10.150.2 5、数字滤波器阶数选择 在数字滤波器设计中，往往需要找到一个满足滤波器设计指标的最小阶数。在 MATLAB 信号处理工具箱中，针对不同的滤波器提供了这样的阶数选择函数。Butterworth 滤波器阶数选择函数 buttord：n,Wn=buttord(Wp,Ws,Rp,Rs)：该函数返回满足通带衰减不大于Rp(dB)、阻带衰减大于Rs(dB)要求的巴特沃思数字滤波器的最小阶数和截止频率。Wp 是通带截止频率，Ws 是阻带截止频率。Wp 和 Ws 的取值都在0 1之间，1 对应的是归一化内奎斯特频率。巴特沃思滤波器的类型有输入参数决定：当 Wp 和 Ws 都是标量并且 WpWs 时，指定滤波器为高通滤波器。通带为（Wp，1），阻带为（0，Ws）；当 Wp 和 Ws 都是二元向量并且Ws(1)Wp(1)Wp(2)Ws(2)时，指定滤波器为带通滤波器，通带为（Wp(1)，Wp(2)），阻带为（0，Ws(1)）和（Ws(2)，1）；当 Wp 和 Ws 都是二元向量并且 Wp(1)Ws(1)Ws(2)Wp(2)时，指定滤波器为带阻滤波器，通带为（0，Wp(1)）和（Wp(2)，1），阻带为（Ws(1)，Ws(2)）。n,Wn=buttord(Wp,Ws,Rp,Rs,s)：这是模拟巴特沃思滤波器的技术选择函数。Chebyshev型滤波器阶数选择函数 cheby1ord 的调用格式如下：n,Wn=cheb1ord(Wp,Ws,Rp,Rs)；n,Wn=cheb1ord(Wp,Ws,Rp,Rs,s)。Chebyshev型滤波器阶数选择函数 cheby2ord 的调用格式如下：n,Wn=cheb2ord(Wp,Ws,Rp,Rs)；n,Wn=cheb2ord(Wp,Ws,Rp,Rs,s)。椭圆滤波器阶数选择函数 ellipord 的调用格式如下：n,Wn=ellipord(Wp,Ws,Rp,Rs)；n,Wn=ellipord(Wp,Ws,Rp,Rs,s)。E.g.设计一带通滤波器，并画出滤波器频率响应。要求，通带波纹 Rp=40dB，通带截止频率 Wp=60 200Hz，阻带截止频率 Ws=50 250Hz，采样频率Fs=1000Hz。（1）MATLAB 程序实现如下：Rp=3;Rs=40;Wp=60 200/500;Ws=50 250/500;N,Wn=cheb2ord(Wp,Ws,Rp,Rs);%切比雪夫型数字滤波器阶数选择函数 N%显示阶数选择值 N=7 Wn%显示截至频率选择值 Wn=0.1000 0.5000 b,a=cheby2(N,Rs,Wn);%切比雪夫型数字滤波器函数 freqz(b,a,512,2000);%绘制频率响应曲线 title(N=,num2str(N);%显示选择阶数（2）程序运行结果 01002003004005006007008009001000-400-2000200400Frequency(Hz)Phase(degrees)01002003004005006007008009001000-