信号处理工具箱

1、1、 信号的表示在信号处理中大多数信号是需要事先输入时间向量的，对于只有一个输入和一个输出的信号，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

2、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

3、0.9894 0.4121 -0.5440使用绘图命令plot(t,y)，可以看到由向量t和y表示的正弦信号如图1.1所示。可见，其图像有些失真，这是因为plot命令是将图中间隔两点用直线连接的，若减小时间间隔p，将有效的恢复正弦函数原貌。>> t=-10:0.01:10;>> y=sin(t);>> plot(t,y)2、 信号的生成1、 正弦波和余弦波在MATLAB中利用函数sin和cos可以生成所需要的正弦波或余弦波。e.g. 生成一个信号持续时长10s，频率为250Hz，幅度为0.75，初始相位为40°的余弦波，并画出其波形图。（1） 问题

4、分析：根据采样定理（采样速率大于等于模拟信号的最高频率的2倍，模拟信号可以由采样序列构成的时间离散信号无失真的表达），采样率至少为信号最高频率的两倍。采样频率必须大于500次/秒，为了产生光滑的曲线，本例中取采样频率为3000次/秒，信号持续时间为10秒，那么采样点数为10×3000=30000。（2） MATLAB命令生成。>> Fs=3000; %采样频率>> t=1/Fs:1/Fs:10; %信号的持续时间>> f=250; %余弦波频率>> A=0.75; %信号幅度>> Ip=40/180*pi; %初始相位>

5、;> y=A*cos(2*pi*f.*t+Ip); %余弦波计算>> plot(t(1:100),y(1:100) %画出余弦波前100个采样值（3） 由余弦波前100点采样值绘出的图形2、 周期方波和锯齿波square命令生成方波，sawtooth命令生成三角波，也称锯齿波。它们的调用格式如下：square(T)：产生一周期为2，幅值为1的方波，采样频率由向量T指定；square(T,DUTY)：产生一个给定占空比，周期为2，幅值为1的方波，占空比是1100之间的数，如果占空比是30，表示一个方波的周期内正电平占30%；sawtooth(T)：产生周期为2，幅值为1的三角波

6、，采样时刻由向量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等份。

7、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');（

8、2） 程序运行结果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的值

9、为负值。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,

10、0.6);>> plot(t,yi,t,yq,t,ye)（2） 程序运行结果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（二次调频）、logarit

11、hmic（对数调频），默认时为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

12、:3; %信号光谱图的实现>> y=chirp(t1,0,1,100);>> specgram(y,256,1e3,256,250); %线性调频信号的光谱图6、 单位冲激信号单位冲激信号是信号系统的基本信号，表示符号(t)，数学定义是： 单位冲激信号除了原点之外，其他各处都为零，并且信号的总面积为1，实际就用一个矩形脉冲来代替单位冲激信号。单位冲激信号在MATLAB中的实现代码如下：>> dt=0.01;>> t=-3:dt:3; %画出-3 3上的波形图>> n=length(t);%计算采样点数n>> x=zeros

13、(1,n);%生成一维数组，其n个元素都为0>> x(1,3/dt+1)=1/dt;%设置原点处的采样值>> stairs(t,x);%画出信号图>> axis(-3,3,0,150)%设置显示窗口横纵坐标 将坐标设置改为 axis(-0.03,0.03,0,150)后可以看到单位冲激信号为窄矩形脉冲近似表示。 7、 单位序列单位序列是一个典型的离散信号，数学表达 。只有当k=0时，函数值才为1，其余全为零，在MATLAB中实现较为简单，下面将给出平移信号(k+k0)的MATLAB实现代码，改信号是指在k=-k0时函数值为1，其余时刻函数值为0 。 e.g

14、画出单位序列 (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(

15、'单位序列')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.

16、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.

17、9998>> x=-4:0.1:6; %画出A中随机数的发生频率分布>> hist(A,x) 9、 高斯分布随机序列在MATLAB中，除了均匀分布的随机序列外，常用的还有标准正态分布的随机序列，该序列可以用randn函数生成，randn函数的调用格式：Y=randn(M,N)，将生成一个M行N列的均值方差为1的标准正态分布的随机数序列。e.g. 产生一个10000个均值为100、方差为4.6的正态分布的随机数序列，并画出其随机数发生频率分布图。MATLAB中实现命令如下：>> M=100;%设置均值>> D=4.6;%设置方差>> Y

18、=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中随机数发生频率分布图3、 随机信号处理和谱分析1、 随机信号互相关函数估计xcorr函数是随机信号互相关估计函数，其调用格式如下。c=xcorr(x,y,maxlags,option)：上式中，x，y为两个独立的随机信号序列，长度都为N的向量（N>1），如果x，y长度不一致，则

19、短的一个用0补齐，使得两个信号长度一样；c为x，y的互相关函数估计序列；maxlags为x与y之间的最大延迟，函数返回值c的长度为2maxlags+1。默认状态下，函数返回值c的长度为2N-1。option指定互相关的归一化选项，它可以是：biased：计算互相关函数的有偏互相关估计；unbiased：计算互相关函数的无偏互相关估计；coeff：系列归一化，使零延迟的自相关为1；none：默认状态，函数执行非归一化计算相关。e.g. 已知两个信号x(t)=sin(2ft)，y(t)=kcos(2ft+)式中，f为10Hz，k为2，为/6。求这两个信号的自相关函数Rx()、Ry()以及互相关函数

20、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(

21、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)的自相

22、关函数>> 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('互相关函数')

23、;>> grid;（2）程序运行结果2、 互协方差函数估计xcov函数是互协方差估计函数，其调用格式如下。c,lags=xcov(x,y,maxlags,option)：上式中，x，y为两个独立的随机信号序列，长度都为N的向量（N>1）；v为x，y的互协方差序列；maxlags为x与y之间的最大延迟，函数返回值c的长度为2N-1。option指定互协方差的归一化选项，它可以是：biased：计算互协方差函数的有偏互相关估计；unbiased：计算互协方差函数的无偏互相关估计；coeff：系列归一化，使零延迟的自相关为1；none：默认状态，函数执行非归一化计算相关。e.g.

24、估计一个正态分布白噪声信号x的自协方差cx(n)，假设最大延迟为30。（1） MATLAB程序实现如下：>> x=randn(1,1000);%生成白噪声>> cov_x,lags=xcov(x,30,'coeff');%计算协方差>> stem(lags,cov_x)%图像输出（2） 程序运行结果3、 谱分析函数psd对于随机信号来说，虽然它没有确定的解析表达式，但其相关函数却是确定的，如果信号是平稳的，那么对相关函数的傅里叶变换就是它的功率谱密度函数，及功率谱，它反映了单位频带内随机信号功率的大小。MATLAB中提供的谱分析函数最常用的是

25、函数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,

26、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；WIND

27、OW为Hanning，数值表述为NFFT；Fs为2；P为0.95；DFLAG为none。其他参数同上。e.g. 采用采样频率为1000Hz、长度为1024点、相邻两段重叠点数为128，窗函数为默认值的Welch方法对信号x(t)=sin2f1t+ksin2f2t+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(

28、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） 程序运行结果4、 谱分析函数pwelchpwelch函数的调用格式如下。Pxx,w=pwelch(x)：该函数用Wel

29、ch方法估计一个输入信号向量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=pwelc

30、h(x,window)：如果设定window是一个正整数，那么这个数代表Hamming窗的长度；如果设定window为一个向量，那么这个向量代表窗函数的权系数。在这种调用格式中，输入向量x被分割成每段重叠50%的整数段，每段的长度和窗函数的长度相同，没有包含在任何一段中的剩余的x中的数据将被忽略。如果指定window为一个向量，则信号数据被分割成8段，并在每一段上加Hamming窗。Pxx,w=pwelch(x,window,noverlap)：改调用格式指定x分割后每一段的长度为window，noverlap指定每段重叠的信号点数，noverlap必须小于被确定的窗口长度，在默认情况下，x被

31、分割后的每段有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为采样频

32、率，如果定义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。twoside

33、d计算双边PSD；onesided计算单边PSD。pwelch(.)：该命令在当前Figure窗口中绘制出功率谱密度曲线，其单位为dB/Hz。e.g. 估计信号x(t)=sin2f1t+ksin2f2t+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;>>

34、; 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） 运行程序结果4、 模拟滤波器设计对特定频率的频点或该频点以外的频率进行有效滤除的系统，就是滤波器，换句话说，滤波器可以被用来区分不同频率的信号，实现各种模拟信号的处理过程。介绍内容：模拟滤波器的设计参数、常用的两种模拟滤波器（巴特沃斯滤波器和切比雪夫滤波器）、模拟滤波器的相互转换。1、 滤波器的设计参数滤波器有四个常用的参数，它们分别是：Fp：通带截

35、止频率，单位Hz；Fs：阻带起始频率，单位Hz；Rs：阻带内最小衰减，又称阻带波纹，单位dB；Rp：通带内允许的最大衰减，又称通带波纹，单位dB。由此，假设采样频率为fN，那么将得到另外两个归一化角频率参数：p:通带截止角频率，单位rad/s，计算公式为p=fp/(fN/2)；s:阻带起始角频率，单位rad/s，计算公式为s=fs/(fN/2)。通过这些参数就可以衡量一个滤波器的性能了，当给定这些参数时，滤波器的设计参数就给定了，用户依照这些参数进行设计。2、 巴特沃斯滤波器巴特沃斯模拟低通滤波器原型的平方幅频响应函数为，c是低通滤波器的截止频率，N为滤波器的阶数。N越大，通带和阻带的近似性越

36、好，过渡带也越陡，特性越接近于矩形。巴特沃思滤波器有以下特点： =c时，A(c2)/A(0)=1/2，幅度衰减，相当于3dB衰减点； 在/c<1的通带内，A（2）有最大平坦的幅度特性，相应（/c）2N随N的增加而趋于0，A（2）趋于1；当/c>1，即在过渡带和阻带中，A（2）单调减小，因为/c>>1，所以A（2）快速下降。在MATLAB中butterworth模拟滤波器函数调用格式为Z,P,K=buttap(N)：N为butterworth滤波器的阶数，Z，P，K分别为滤波器的零点、极点、增益。该函数返回N阶低通模拟滤波器原型的极点和增益。e.g. 画出5阶和10阶bu

37、tterworth低通滤波器的平方幅频响应曲线。（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);

38、>> 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） 程序运行结果3、 切比雪夫型滤波器巴特沃斯滤波器的频率特性无论在通带与阻带都随频率而单调变化，因而如果在通带边缘满足指标，则在通带内肯定满足，也就是会超过指标的要求，因而并不经济。所以，更有效的办法是将指标的精度要求均匀地分布在通带内，或均匀

39、分布在阻带内，或同时均匀地分布在通带与阻带内，这样在同样通带、阻带性能要求下，就可设计出阶数较低的滤波器。这种滤波器均匀分布的办法可通过选择具有等波纹特性的逼近函数来完成。切比雪夫滤波器的幅度特性就在一个频带中（通带或阻带）具有这种等波纹特性。切比雪夫型滤波器的平方幅频响应函数为：。式中，为小于1的正数，表示通带内幅频波纹情况；c是截止频率，N为多项式的阶数，其中。MATLAB中切比雪夫型滤波器函数为：Z,P,K=cheb1ap(N,Rp)：N为阶数，Z、P、K分别为滤波器的零点、极点、增益；Rp为通带波纹。e.g. 画出5阶切比雪夫型模拟低通滤波器原型的平方幅频响应曲线。（1） MATLAB

40、程序实现如下：>> 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')>>

41、; title('N=5');（2） 程序运行结果4、 切比雪夫型滤波器切比雪夫型滤波器的平方幅频响应函数为。式中，为小于1的正数，表示阻带内幅频波纹情况；c是截止频率，N为多项式的阶数，其中。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=che

42、b2ap(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） 程序运行结果5、 模拟滤波器的

43、频域变换如果想设计其他类型的频率选择性滤波器，如高通、带通、带阻滤波器，可以将一个低通滤波器的频带进行变换，以使其能够表现其他频率选择性滤波器特性。 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为

44、转换后高通滤波器的截止频率，函数将模拟低通滤波器转换为指定带宽和中心频率的高通滤波器。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)：将连续时间状态空间低通滤波器原型转换

45、为一个中心频率为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=32

46、0Hz、带宽为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);%模拟低通滤波器转换为低通滤波器>&

47、gt; 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;>> fre

48、qs(num2,den2);%绘制频率响应曲线5、 IIR数字滤波器设计数字滤波器是数字信号处理的基础，用来对信号进行过滤、检测与参数估计等处理，在通信、图像、语音、雷达等许多领域都有着十分广泛的应用目前数字滤波器的设计有许多现成的高级语言设计程序，但它们都存在设计效率较低，不具有可视图形，不便于修改参数等缺点，而MATLAB为数字滤波的研究和应用提供了一个直观、高效、便捷的通道。MATLAB中自带了许多IIR数字滤波器设计函数。1、 巴特沃思数字滤波器设计（butter函数）butter函数即可以用来设计低通、带通、高通和带阻的模拟巴特沃斯滤波器，也可以用来设计这四类的数字巴特沃思滤波器。数

49、字滤波器设计的butter函数调用格式如下： b,a=butter(n,Wn)：函数用来设计一个截止频率为Wn的n阶低通滤波器。它返回滤波器系数向量a、b的长度为n+1，这些系数按z的降幂排列为：截止频率是滤波器的幅度响应为时的频率，归一化截止频率Wn取值在0 1之间，这里1对应内奎斯特频率。如果Wn是一个二元向量，如Wn=w1 w2；那么函数返回一个带通为w1<w<w2、阶数为2×n的带通数字滤波器。b,a=butter(n,Wn,ftype)：该函数格式设计一个截止频率为Wn的高通或者带通数字滤波器，ftype为滤波器类型参数，high为高通滤波器；stop为带阻滤波

50、器。对于带阻滤波器，如果Wn是一个二元向量，如Wn=w1 w2，那么butter函数返回一个带通为w1<w<w2、阶数为2×n的带通数字滤波器。z,p,k=butter(n,Wn)或z,p,k=butter(n,Wn,ftype)：这是butter函数的零极点形式，其返回零点和极点的n列向量z、p，以及增益标量k，其他参数同上。A,B,C,D=butter(n,Wn)或A,B,C,D=butter(n,Wn,ftype)：这是butter函数的状态空间形式，这里的A、B、C、D的关系如下：x(n+1)=Ax(n)+Bu(n)；y(n)=Cx(n)=Du(n)这里u是输入，

51、y是输出，x是状态向量，其他参数同上。模拟滤波器设计的butter函数调用格式为：b,a=butter(n,Wn,'s')b,a=butter(n,Wn,'ftype','s')z,p,k=butter(n,Wn,'s')z,p,k=butter(n,Wn,'ftype','s')A,B,C,D=butter(n,Wn,'s')A,B,C,D=butter(n,Wn,'ftype','s')e.g. 设计一个9阶的巴特沃思高通滤波器，采样频率为1000

52、Hz，截止频率为300Hz（归一化频率值为0.6），并画出滤波器频率响应曲线。（1） MATLAB程序实现如下>> N=9;>> Wn=300/500;>> b,a=butter(N,Wn,'high');%巴特沃思高通数字滤波器函数>> freqz(b,a,128,1000);%绘制频率响应曲线>> title('N=9');%显示阶数（2） 程序运行结果e.g. 设计一个10阶的巴特沃思带通滤波器，通带频率为100Hz到200Hz，并画出其脉冲响应曲线。（1） MATLAB程序实现如下>>

53、; n=5;>> Wn=100 200/500;>> b,a=butter(n,Wn);>> y,t=impz(b,a,101);>> stem(t,y)（2） 程序运行结果2、 切比雪夫型数字滤波器设计（cheby1函数）cheby1函数即可以用来设计低通、带通、高通和带阻的模拟切比雪夫型滤波器，也可以用来设计这四类的数字切比雪夫型滤波器。它的特性就是通带内等波纹，阻带内单调。cheby1函数调用格式如下。 b,a=cheby1(n,Rp,Wn)：函数用来设计一个截止频率为Wn，通带波纹为Rp(dB)的n阶切比雪夫低通滤波器。它返回滤波器系数向

54、量a、b的长度为n+1，这些系数按z的降幂排列为：归一化截止频率是滤波器的幅度响应为-Rp(dB)时的频率，对于cheby1来说，归一化截止频率Wn取值在0 1之间，这里1对应内奎斯特频率。如果Wn是一个二元向量，如Wn=w1 w2，那么cheby1函数返回一个带通为w1<w<w2、阶数为2×n的带通数字滤波器。b,a=cheby1(n,Rp,Wn,ftype)：该函数格式设计一个截止频率为Wn、通带波纹为Rp(dB)的高通或者带通数字滤波器，ftype为滤波器类型参数，high为高通滤波器；stop为带阻滤波器。对于带阻滤波器，如果Wn是一个二元向量，如Wn=w1,w2

55、，那么cheby1函数返回一个带通为w1<w<w2、阶数为2×n的带通数字滤波器。z,p,k=cheby1(n,Rp,Wn)或z,p,k=cheby1(n,Rp,Wn,ftype)：这是cheby1函数的零极点形式，其返回零点和极点的n列向量z、p，以及增益标量k。其他参数同上。A,B,C,D=cheby1(n,Rp,Wn)或A,B,C,D=cheby1(n,Rp,Wn,ftype)：这是cheby1函数的状态空间形式，这里的A、B、C、D的关系如下：xn+1=Ax(n)+Bu(n)；y(n)=Cx(n)+Du(n)这里u是输入，y是输出，x是状态向量。其他参数同上。模拟

56、滤波器设计的cheby1函数调用格式为：b,a=cheby1(n,Rp,Wn,'s')b,a=cheby1(n,Rp,Wn,'ftype','s')z,p,k=cheby1(n,Rp,Wn,'s')z,p,k=cheby1(n,Rp,Wn,'ftype','s')A,B,C,D=cheby1(n,Rp,Wn,'s')A,B,C,D=cheby1(n,Rp,Wn,'ftype','s')e.g. 设计一个9阶的cheby1型低通数字滤波器，采样频率为10

57、00Hz，Rp=0.5dB，截止频率为300Hz（归一化频率值为0.6），并画出滤波器频率响应曲线。（1） MATLAB程序实现如下：>> 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） 程序运行结果e.g. 设计一个10阶的切比雪夫型带通滤波器，通带频率为

58、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） 程序运行结果3、 切比雪夫型数字滤波器设计（cheby2函数）Cheby2函数即可以用来设计低通、带通、高通和带阻的模拟切比雪夫型滤波器，也可以用来设计这四类的数字切比雪夫型滤波器。它的特性就是通带内单调，阻带内等波纹。数字滤波器设计的ch

59、eby2函数调用格式为：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,

60、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） 程序运行结果e.g. 设计一个10阶的切比雪夫型带通滤波器，