《数字信号处理》上机实验指导

《数字信号处理》上机实验指导 实验一、Z变换及离散时间系统 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 （一）、实验目的 1、通过本实验熟悉Z变换在离散时间系统分析中的地位和作用。 2、掌握并熟练使用有关离散系统分析的MATLAB调用函数及格式，以深入理解离散时间系统的频率特性。 （二）、实验内容及步骤 对于一个给定的LSI系统，其转移函数H(z)习惯被定义为H(z)=B(z)/A(z),即： 公式中 和 分别是H(Z)分子与分母多项式的阶次，在有关MATLAB 的系统分析的文件中，分子和分母的系数被定义为向量，即 并要求 =1,如果 ≠1，则程序将自动的将其归一化为1。 1、系统的阶跃响应 调用格式为： y=filter(b,a,x),其中x,y,a,b都是向量。 例1 令 求该系统的阶跃响应（y（n））。实现该任务的程序如下： clear; x=ones(100);% x(n)=1,n=1~100； t=1:100;% t用于后面的绘图； b=[.001836，.007344，.011016，.007374，.001836]; % 形成向量b； a=[1，-3.0544，3.8291，-2.2925，.55075]; % 形成向量a； y=filter(b,a,x); % 求所给系统的输出，本例实际上是求所给系统的阶跃响应； plot(t,x,'r.',t,y,'k-');grid on;% 将x(n)（绿色）y(n)（黑色）画在同一个 %图上； ylabel('x(n) and y(n)') xlabel('n') 2、单位抽样响应h(n) 调用格式为： h=impz(b，a，N) 或 [h，t]=impz(b，a，N) 其中N是所需的h(n)的长度，前者绘图时n从1开始，而后者从0开始。 例2、求上例所给系统的单位抽样响应h(n)。 实现该任务的程序如下： clear; b=[.001836，.007344，.011016，.007374，.001836]; a=[1，-3.0544，3.8291，-2.2925，.55075]; [h,t]=impz(b，a，40); % 求单位抽样响应 stem(t，h，'.');grid on; 3、求频率响应 基本调用格式为：[H，w]=freqz(b，a，N,‘whole’,Fs) 其中N是频率轴的分点数，建议N为2的整次幂；w是返回频率轴坐标向量，供绘图用；Fs是抽样频率，若Fs=1,频率轴给出归一化频率；whole指定计算的频率范围是从0~Fs,缺省时是从0~Fs/2。 例3、求例1所给系统的频率响应 。 实现该任务的程序如下： clear all; b=[.001836,.007344,.011016,.007374,.001836]; a=[1,-3.0544,3.8291,-2.2925,.55075]; [H,w]=freqz(b,a,256,1); Hr=abs(H);%绝对值（幅值）； Hphase=angle(H);%相位角； Hphase=unwrap(Hphase); % 解卷绕 subplot(211) plot(w,Hr);grid on; ylabel(' Amplitude Freq. Res.') subplot(212) plot(w,Hphase);grid on; ylabel(' Phase Freq. Res.') 4、离散系统的极零图 调用格式：zplane(z,p) 或zplane(b,a) 前者是在已知系统零点的列向量z和极点的列向量p的情况下画出的极零图，后者是在已知B(z)，A(z)的情况下的极零图。 例4、显示例1系统及FIR系统 的极零图。 实现该任务的程序为： clear; b=[.001836,.007344,.011016,.007374,.001836]; a=[1,-3.0544,3.8291,-2.2925,.55075]; subplot(221); zplane(b,a); % 求并画出所给系统的极零图，该系统为IIR系统； b=[1 -1.7 1.53 -0.68]; a=1; subplot(222); zplane(b,a); % 求并画出第二个系统的极零图，该系统为FIR系统 （三）、作业 给定系统 ，编程并绘出系统的单位阶跃响应y(n),频率响应 。给出实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 。 clear; x=ones(100);% x(n)=1,n=1~100； t=1:100;% t用于后面的绘图； b=[0,0,-0.2]; % 形成向量b； a=[1,0,0.8]; % 形成向量a； y=filter(b,a,x); % 求所给系统的输出，本例实际上是求所给系统的阶跃响应； subplot(311) plot(t,x,'r.',t,y,'k-');grid on;% 将x(n)（绿色）y(n)（黑色）画在同一个 %图上； ylabel('x(n) and y(n)') xlabel('n') [H,w]=freqz(b,a,256,1); Hr=abs(H);%绝对值（幅值）； Hphase=angle(H);%相位角； Hphase=unwrap(Hphase); % 解卷绕 subplot(312) plot(w,Hr);grid on; ylabel(' Amplitude Freq. Res.') subplot(313) plot(w,Hphase);grid on; ylabel(' Phase Freq. Res.') 实验二、快速傅里叶变换 （一）、实验目的 1、通过本实验进一步加深对快速傅里叶变换的理解。 2、会熟练运用fft,ifft,czt实现线性调频z变换。 （二）、实验内容 1、快速傅里叶变换（fft） 调用格式为 X=fft(x) 或 X=fft(x,N) 对前者，若x的长度是2的整次幂，则按该长度实现x的快速变换，否则，实现的是非2 的整次幂的变换；对后者，N应为2的整次幂，若x得长度小于N，则补零，若超过N，则舍弃N以后的数据。ifft的调用格式与之相同。 例1、令x(n)是两个正弦信号及白噪声的叠加，试用fft文件对其作频谱分析。 实现该任务的程序为： clear all; % 产生两个正弦加白噪声； N=256;% 2的8次幂，进行8级蝶形运算； f1=.1;f2=.2;fs=1; a1=5;a2=3; w=2*pi/fs; x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N); % 应用FFT 求频谱； subplot(3,1,1); plot(x(1:N/4));%作图向量x的第一到第N/4个值； f=-0.5:1/N:0.5-1/N; X=fft(x); % 快速傅里叶变换； y=ifft(X);% 快速傅里叶逆变化 subplot(3,1,2); plot(f,fftshift(abs(X))); subplot(3,1,3); plot(real(y(1:N/4))); %作图向量y的第一到第N/4个值的实部； 附:fftshift函数 fftshift函数就是一个交换函数： 交换 规则 编码规则下载淘宝规则下载天猫规则下载麻将竞赛规则pdf麻将竞赛规则pdf 如下： 如:x=[1 2 3 4 5 6 7 8]; y=fftshift(x); then y=[5 6 7 8 1 2 3 4]; 其在fft运算里的物理意义： 把0频(低频)周围的频谱搬移到中频范围(采样频率的一半)，只是形象化的展示FT变换后的低频成分(正负频率)。 其实质是把Fs/2的右边频谱平移到Fs/2的左边，把低频平移到Fs/2的右边，各图象间距不变。 2、线性调频Z变换（CZT） CZT可用来计算单位圆上任一段曲线上的Z变换，做DFT时输入的点数N和输出的点数可以不相等，从而达到频域“细化”的目的。CZT在单位圆上的Z变换就是傅里叶变换。 其调用格式为： X=czt(x,M,W,A) 式中x是待变换的时域信号x(n)，其长度设为N，M 是变换的长度，W确定变换的步长，A确定变换的起点。若M=N,A=1，则CZT变成DFT。CZT应用举例： 例2、设x(n)由三个实正弦组成，频率分别是8Hz,8.22Hz,9Hz,抽样频率为40 Hz,时域取128点，作CZT变换和FFT变换。 实现该任务的程序为： clear all； % 构造三个不同频率的正弦信号的叠加作为试验信号 N=128; f1=8;f2=8.22;f3=9;fs=40; stepf=fs/N; n=0:N-1; t=2*pi*n/fs; n1=0:stepf:fs/2-stepf; x=sin(f1*t)+sin(f2*t)+sin(f3*t); M=N; W=exp(-j*2*pi/M); % A=1时的czt变换 A=1; Y1=czt(x,M,W,A); subplot(311) plot(n1,abs(Y1(1:N/2)));grid on; % DTFT Y2=abs(fft(x)); subplot(312) plot(n1,abs(Y2(1:N/2)));grid on; % 详细构造A后的czt M=60; f0=7.2; DELf=0.05; A=exp(j*2*pi*f0/fs); W=exp(-j*2*pi*DELf/fs); Y3=czt(x,M,W,A); n2=f0:DELf:f0+(M-1)*DELf; subplot(313);plot(n2,abs(Y3));grid on; （三）、作业 设x(n)由三个实正弦组成，频率分别是8Hz,9Hz,10Hz,抽样频率为60 Hz,时域取256点，作CZT变换、IFFT变换和FFT变换，观察波形，更改参数，得出不同参数下的CZT变换波形。给出实验报告。 clear all; N=256;% 2的8次幂，进行8级蝶形运算 f1=3;f2=5;f3=7;fs=40; a1=1;a2=3;a3=5; stepf=fs/N; n=0:N-1; t=2*pi*n/fs; n1=0:stepf:fs/2-stepf; x=a1*sin(f1*t)+a2*sin(f2*t)+a3*sin(f3*t); M=N; W=exp(-j*2*pi/M); % 三个不同频率的正弦信号； subplot(511); plot(n,x);grid on; ylabel('one') % 应用FFT 求频谱； X=fft(x); % 快速傅里叶变换； subplot(512); plot(n1,abs(X(1:N/2)));grid on; ylabel('two') y=ifft(X);% 快速傅里叶逆变化 subplot(513); plot(real(y(1:N)));grid on; ylabel('three') % A=1时的czt变换 A=1; Y1=czt(x,M,W,A); subplot(514) plot(n1,abs(Y1(1:N/2)));grid on; ylabel('four') % 详细构造A后的czt M=60; f0=7.2; DELf=0.05; A=exp(j*2*pi*f0/fs); W=exp(-j*2*pi*DELf/fs); Y3=czt(x,M,W,A); n2=f0:DELf:f0+(M-1)*DELf; subplot(515);plot(n2,abs(Y3));grid on; ylabel('five') 实验三、无限冲击响应数字滤波器设计 （1）​ 实验目的 1、​ 要求掌握IIR数字滤波器的设计原理、设计 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 和设计步骤； 2、​ 能根据给定的滤波器指标进行滤波器设计； 3、​ 掌握数字巴特沃斯滤波器、数字切比雪夫滤波器的设计原理和步骤； （二）、实验内容 IIR数字滤波器的设计有多种方法，如频率变换法、数字域直接设计以及计算辅助设计等。下面只介绍频率变换设计法。首先考虑由模拟低通滤波器到数字低通滤波器的转换，其基本的设计过程如下: 1、​ 将数字滤波器的技术指标转换为模拟滤波器的技术指标； 2、​ 设计模拟滤波器G(S)； 3、​ 将G(S)转换成数字滤波器H(Z) 在低通滤波器的设计基础上，可以得到数字高通、带通、带阻滤波器的设计 流程 快递问题件怎么处理流程河南自建厂房流程下载关于规范招聘需求审批流程制作流程表下载邮件下载流程设计 ，以高通数字滤波器的设计为例： 1、​ 将数字高通滤波器 的技术指标 ，通过 转变为模拟高通 的技术指标 ，作归一化处理后得 ； 2、​ 利用频率变换关系 ，将模拟高通 的技术指标转换为归一化的低通滤波器G(p)的技术指标，并有p= ; 3、​ 设计模拟低通滤波器G(p); 4、​ 将G(p)转换为模拟高通滤波器的转移函数 ,p= ; 5、​ 将 转换成数字高通滤波器的转移函数 ，s=(z-1)/(z+1)。 以上5个步骤同样适用于数字带通、数字带阻滤波器的设计。只是在步骤2，3，4中频率转换的方法不同。 例1、试设计一个数字高通滤波器，要求带通下限频率 =0.8pi,阻带上限频率为0.44pi,通带衰减不大于3dB,阻带衰减不小于20 dB。 实现该任务的MATLAB程序为： clear all wp=.8*pi;ws=.44*pi; rp=3;rs=20;% 通带和阻带的衰减 Fs=2000;% 抽样衰减 %数字高通滤波器的技术指标 % Firstly to finish frequency prewarping; wap=2*Fs*tan(wp/2) was=2*Fs*tan(ws/2);% 模拟高通滤波器的技术指标 [n,wn]=buttord(wap,was,rp,rs,'s'); % 确定模拟低通滤波器的阶数 % Note: 's'! [z,p,k]=buttap(n);% 设计模拟低通原型滤波器，z,p,k分别是设 计出低通原型滤波器的极点，零点及增益 [b,a]=zp2tf(z,p,k)% 零、极点系数转换为传递函数 % 巴特沃斯模拟低通原型滤波器频率响应 [bt,at]=lp2hp(b,a,wap)%由低通滤波器转移函数产生模拟高通滤波器转移函数 % % Note: z=(2/ts)(z-1)/(z+1);2/ts=1,that is 2fs=1,fs=0.5; [bz,az]=bilinear(bt,at,Fs)% 由模拟高通滤波器H(s)得到数字高通滤波器H(z) [h,w]=freqz(bz,az,256,1);%得到频率响应 figure(1) plot(w,20*log10(abs(h)));grid;%得出幅频图 ylabel(' High-pass DF') % % Directly to design H(z) by butter.m [n,wn]=buttord(wp/pi,ws/pi,rp,rs);% 确定模拟低通滤波器的阶数 [b,a]=butter(n,wp/pi,'high');%直接设计巴特沃思滤波器 [h1,w1]=freqz(b,a,256,1); h1=20*log10(abs(h1)); figure(2) plot(w,h1);grid ylabel(' High-pass DF') 例2、一个数字系统的抽样频率为1000Hz,已知该系统受到频率为100Hz的噪声的干扰，现设计一陷波滤波器 去掉该噪声。要求3dB的带边频率为95Hz和105Hz,阻带衰减不小于14Db.阻带的下边和上边频率分别为99Hz和101Hz。 实现程序为： fp=[95 105];fs=[99 101]; %wp=[.19*pi 0.21*pi];ws=[.198*pi 0.202*pi]; Fs=1000; rp=3;rs=14; wp=fp*2*pi/Fs;ws=fs*2*pi/Fs; % % Firstly to finish frequency prewarping; wap=2*Fs*tan(wp./2) was=2*Fs*tan(ws./2); [n,wn]=buttord(wap,was,rp,rs,'s'); % Note: 's'! [z,p,k]=buttap(n); [b,a]=zp2tf(z,p,k) bw=wap(2)-wap(1) w0=sqrt(wap(1)*wap(2)) [bt,at]=lp2bs(b,a,w0,bw)%将模拟低通滤波器G(p)转换为实际带阻滤波器 % % Note: z=(2/ts)(z-1)/(z+1); [bz1,az1]=bilinear(bt,at,Fs) [h,w]=freqz(bz1,az1,256,Fs); figure(1) plot(w,20*log10(abs(h))) % Directly to design H(z) by butter.m [n,wn]=buttord(wp/pi,ws/pi,rp,rs); [b,a]=butter(n,wp/pi,'stop') [h1,w1]=freqz(b,a,256,Fs); figure(2) plot(w1,20*log10(abs(h1))) 附： buttord.m 本文件用来确定数字低通或模拟低通滤波器的阶数，其调用格式为 (1)​ [N,Wn]= buttord(Wp, Ws, Rp, Rs) (2)​ [N,Wn]= buttord(Wp, Ws, Rp, Rs,‘s’) 格式（1）对应数字滤波器，式中Wp, Ws分别是通带和阻带的截止频率，实际上他们是归一化频率，其值在0~1之间，1对应抽样频率的一半。Rp, Rs分别是通带和阻带的衰减，单位是Db。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率。格式（2）对应模拟滤波器，式中各个变量的含义与格式1相同，但Wp, Ws及Wn的单位为rad/s,因此它们实际上是频率 。 butter.m 本文件可用来直接设计巴特沃思数字滤波器。其调用格式为 （1）[B,A]=butter（N，Wn）; （2）[B,A]=butter（N，Wn，‘high’）; （3）[B,A]=butter（N，Wn，‘stop’）; （4）[B,A]=butter（N，Wn，‘s’）; 格式1~3用来设计数字滤波器,所以B,A分别是H(z)的分子，分母的系数向量，Wn是通带截止频率，其值在0~1之间，1对应抽样频率的一半。若Wn是标量，则格式1用来设计低通数字滤波器，若Wn是1*2的向量，则格式1用来设计带通数字滤波器，2用来设计高通数字滤波器，3用来设计带阻数字滤波器，4用来设计模拟滤波器。 lp2bp.m 本文件用来将模拟低通滤波器G(p)转换为实际的带通滤波器，其调用格式为 [B,A]= lp2bp(b,a,WO, BW) 式中b,a分别是模拟低通原型滤波器G(p)的分子、分母多项式的系数向量，B,A分别是转换后的H(s)的分子、分母多项式的系数向量， WO 是带通滤波器的中心频率，BW是其带宽。 （三）、作业 设计一个数字带通滤波器，参数自定。给出实验报告。clear all; fp=[100 400];fs=[50 450]; Fs=1000;% 抽样衰减 rp=3;rs=18;% 通带和阻带的衰减 %数字高通滤波器的技术指标 wp=fp*2*pi/Fs;ws=fs*2*pi/Fs; wap=2*Fs*tan(wp./2) was=2*Fs*tan(ws./2); [n,wn]=buttord(wap,was,rp,rs,'s'); [z,p,k]=buttap(n);% 设计模拟低通原型滤波器，z,p,k分别是设计出低通原型滤波器的极点，零点及增益 [b,a]=zp2tf(z,p,k)% 零、极点系数转换为传递函数 % 巴特沃斯模拟低通原型滤波器频率响应 bw=wap(2)-wap(1) w0=sqrt(wap(1)*wap(2)) [bt,at]=lp2bp(b,a,w0,bw) [bz1,az1]=bilinear(bt,at,Fs) [h,w]=freqz(bz1,az1,256,Fs);%得到频率响应 figure(1) plot(w,20*log10(abs(h))) grid on 实验四、有限冲击响应数字滤波器设计 (一)、实验目的 1、了解无限冲击响应数字滤波器设计和有限冲击响应数字滤波器 设计各自的特点，比较两者的优缺点。 2、掌握用窗函数法设计FIR 数字滤波器的原理和方法。 3、熟悉线性相位FIR 数字滤波器特性。 4、了解各种窗函数对滤波特性的影响。 （二）、实验原理及内容 窗口法的基本思想：根据给定的滤波器技术指标，选择滤波器长度N和窗函数 ，使其具有最窄宽度的主瓣和最小的旁瓣。 1、窗口法的设计步骤： （1）、给定理想频响函数 ； （2）、根据指标选择窗函数。确定窗函数类型的主要依据是过度带宽和阻带最小衰耗的指标； （3）、由 求 ，加窗得 ； （4）、检验。由 求 ，求 是否在误差容限之内。 2、常用MATLAB文件介绍 fir1.m 本文件采用窗函数法设计FIR数字滤波器，其调用格式为 (1) b=fir1(N,Wn) (2) b=fir1(N,Wn,‘high’) (3) b=fir1(N,Wn,‘stop’) 式中N为滤波器的阶数，因此滤波器的长度为N+1; Wn是通带截止频率，其值在0~1之间，1对应抽样频率的一半；b是设计好的滤波器系数h(n).对格式1，若Wn是一标量，则可用来设计低通滤波器；若Wn是1*2的向量，则用来设计带通滤波器；若Wn是1*L的向量，则可用来设计L带通滤波器，这时，格式1 要改为 b=fir1(N,Wn,‘DC-1’) 或 b=fir1(N,Wn,‘DC-0’) 前者保证第一个带为通带，后者保证第一个带为阻带。格式2用来设计高通滤波器，3用来设计带阻滤波器。在上述所有格式中，若不指定窗函数的类型，fir1自动选择汉明窗。 例1、​ 设计一FIR低通滤波器，所希望的频率响应 在0≤ ≤0.25pi之间，在0.25pi≤ ≤pi之间为0，令N=10，分别用矩形窗和汉明窗设计，观察其幅频响应的特性。 实现该任务的程序为： clear all; N=10; b1=fir1(N,0.25,boxcar(N+1)); % 用矩形窗作为冲激响应的窗函数 b2=fir1(N,0.25,hamming(N+1));% 用Hamming窗作为冲激响应的窗函数 M=128; h1=freqz(b1,1,M); h2=freqz(b2,1,M); % 分别求两个滤波器的频率响应； t=0:10; subplot(221) stem(t,b2,'.');hold on;% 绘制火柴梗图； plot(t,zeros(1,11));grid;%绘制1*11的零数组网格图； f=0:0.5/M:0.5-0.5/M; subplot(222) plot(f,abs(h1),'b-',f,abs(h2),'g-');grid; fir2.m 本文件采用窗函数法设计具有任意幅频响应的FIR数字滤波器，其调用格式为 b=fir2(N,F,M) 其中F是频率向量，其值在0~1之间，M是和F相对应的所希望的幅频响应。如同fir1缺省时自动选用汉明窗。 例2、​ 设计一多带滤波器，要求理想幅频响应在归一化频率 0.2~0.3，0.6~0.8之间为1，其余处为0。 实现程序如下（给出了两个不同长度滤波器以作比较）： clear all; f=[0 0.19 0.2 0.3 0.31 0.59 0.6 0.8 0.81 1]; % 给定频率轴分点； m=[0 0 1 1 0 0 1 1 0 0]; % 给定在这些频率分点上理想的幅频响应 N1=30; N2=90; % 取两种不同的滤波器长度； b1=fir2(N1,f,m); b2=fir2(N2,f,m); % 得到两个滤波器； subplot(311); stem(b1,'.');grid; subplot(312); stem(b2,'.');grid; M=128; [h1,w]=freqz(b1,1,M,1); [h2,w]=freqz(b2,1,M,1); subplot(313); plot(w,abs(h1),'b-',w,abs(h2),'g-');grid; （三）、作业 设计一FIR低通滤波器，所希望的频率响应 在0≤ ≤0.3pi之间，在0.3pi≤ ≤pi之间为0，分别取N=10，20，40，自行选择窗函数，观察其幅频响应的特性。给出实验报告。 clear all; N=10; b1=fir1(N,0.3,boxcar(N+1)); b2=fir1(N,0.3,hamming(N+1)); M=128; h1=freqz(b1,1,M); h2=freqz(b2,1,M); t=0:10; subplot(221) stem(t,b2,'.');hold on; plot(t,zeros(1,11));grid; f=0:0.5/M:0.5-0.5/M; subplot(222) plot(f,abs(h1),'b-',f,abs(h2),'g-');grid; clear all; N=20; b1=fir1(N,0.3,boxcar(N+1)); b2=fir1(N,0.3,hamming(N+1)); M=128; h1=freqz(b1,1,M); h2=freqz(b2,1,M); t=0:20; subplot(221) stem(t,b2,'.');hold on; plot(t,zeros(1,21));grid; f=0:0.5/M:0.5-0.5/M; subplot(222) plot(f,abs(h1),'b-',f,abs(h2),'g-');grid; clear all; N=40; b1=fir1(N,0.3,boxcar(N+1)); b2=fir1(N,0.3,hamming(N+1)); M=128; h1=freqz(b1,1,M); h2=freqz(b2,1,M); t=0:40; subplot(221) stem(t,b2,'.');hold on; plot(t,zeros(1,41));grid; f=0:0.5/M:0.5-0.5/M; subplot(222) plot(f,abs(h1),'b-',f,abs(h2),'g-');grid;