数字信号处理报告(20211125224422)

1、数字信号处理实验报告实验目的一、加深对离散傅立叶变换（dft）和快速傅立叶变换（fft）的理解，掌握通过两种变换求卷积的编程方法；二、掌握设计巴特沃斯低通双线性iir 数字滤波器的原理和方法，以及从低通转换到高通的技术；三、掌握用窗函数设计fir 数字滤波器的原理及方法，了解各种不同窗函数对滤波器性能的影响；四、提高综合应用和分析的能力，matlab 编程能力等。实验内容实验一一、实验名称快速傅立叶变换二、实验要求编程利用fft 进行卷积计算，通过实验比较出快速卷积优越性。三、实验原理利用 fft 进行离散卷积的步骤归纳如下：(1)、设 x(n)的列长为 n1，h(n)的列长为 n2，要求y(

2、n)=x(n)*h(n)=10nkx(k)h(n-k)nrn(n)=10)()(nkknhkx1 (2)、为使两有限长序列的线性卷积可用其圆周卷积来代替而不产生混淆，必须选择nn1+ n2-1 。为使用基 -2fft 来完成卷积计算，故要求n=2v(v 是整数 ) 。用补零的办法使x(n) ，h(n) 具有列长 n，即x(n)=1-n1,n1n1n01-1210nx(n)，nh(n)=1-n1,n2n2n01-2210nh(n)，n(3)为用圆周卷积定理计算线性卷积，先用fft 计算 x(n) ，h(n) 的 n 点离散傅立叶变换x(n)fftx(k) 2 h(n) ffth(k) 3 (4)

3、 组成卷积y(k)=x(k)h(k) 4 (5)利用 ifft 计算 y(k) 的离散傅立叶逆变换得到线性卷积y(n)。由于y(n)=10)()/1(nkkynwn-nk=10)(*)/1(nkkynwnnk* 5 可见， y(n)可由求 (1/n)y*(k) 的 fft 再取共轭得到。四、实验题目(1)两个正弦序列的卷积（均为两个周期，256 点）输入序列：卷积输出：（2）正弦序列与三角序列的卷积（正弦序列为两个周期，256 点；三角序列为一个周期，256 个点）输入序列：卷积输出：（3）两个矩形序列的卷积（均为两个周期，256 个点，占空比0.5）输入序列：卷积输出：（4）单位冲击与正弦波

4、（单位冲击序列为256 个点，正弦序列为1.7 个周期， 256 个点）输入序列：卷积输出：任何序列与单位冲击序列的卷积为原序列，所以结果正确。(5)正弦序列与矩形序列的卷积（两序列均为256 个点，正弦序列为两个周期，矩形序列为两个周期，占空比为0.2）输入序列：卷积输出：主要代码（只选一个参考，下同）：% （1）的代码k=1:256; s1=sin(k/64*pi); s2=s1; xk=fft(s1,2*length(k)-1); yk=fft(s2,2*length(k)-1); rm=ifft(xk.*yk); m=(-255):(255); stem(m,rm) xlabel(m)

5、; ylabel( ? ); 实验二一、实验名称iir 数字滤波器二、实验要求：设计巴特沃斯低通双线性iir 数字滤波器， n=3，c=0.2 。输入信号为以下三种序列（选择点数为128 或 256） ：1 单位取样2 三角序列（两个周期）3 矩形序列（占空比0.1 或 0.5 ）给出输出的时域及频域效果，并进行简单的分析。三、实验原理：滤波器的作用是滤除信号中某一部分的频率分量。信号经过滤波器处理，就相当与信号频谱与滤波器的频率响应相乘的结果。在时域里来看，这就是信号与滤波器的冲激响应相卷积。可以说滤波器就是一个卷积器。iir 滤波器的系统函数对应的差分方程为模拟滤波器系统函数)(sha的一

6、般表示式为mrnkkrknyarnxbny01)()()(nkkkmrrrzazbzh101)(数字滤波器系统函数h(z) 的普遍表示式为三阶次巴特沃斯滤波器的系统函数为由 ha(s) 的系数表示经双极性变换后的y(z) 的表达式（三阶）c=0.2 y(z)=0.018099*x(z)*z-3+0.054297*x(z)*z-2+0.054297*x(z)*z-1+0.018099*x(z) +0.27806*y(z)*z-3-1.18289*y(z)*z-2+1.76004*y(z)*z-1即y(n)=0.018099*x(n-3)+0.054297*x(n-2)+0.054297*x(n-

7、1)+0.018099*x(n) +0.27806*y(n-3)-1.18289*y(n-2)+1.76004*y(n-1) 由低通数字滤波器原形变换为高通数字滤波器由截止频率 c=0.2 三阶低通变换为截止频率c=0.6 三阶高通，经计算，其表达式为：y(n)=-0.098531x(n-3)+0.295594x(n-2)-0.295594x(n-1)+0.098531x(n) -0.056297y(n-3)-0.42179y(n-2)-0.57724y(n-1) 计算过程：由给定的条件可计算出巴特沃斯系统函数的系数，相应可知摸拟系统函数的系数，经双极性变换法求出数字滤波器的系数，最后由差分方

8、程实现低通滤波效果。经相应的z 平面映射，由映射公式变换得出数字高通滤波器系统函数的系数，从而由差分方程实现高通效果。四、实验结果及分析：a.低通滤波：)22/(32233ccccsss)2tan(2tcnnnnnkkknrrrascscsccsdsdsddscsdsh2210221000)(nnnnnkkknrrrzazazazbzbzbbzazbzh221122110001)(1111aaz)2cos()2cos(cccca(1). 单位取样（ 256 个点）输入序列：低通滤波后的傅立叶变换和输出序列：0100200300-0.500.51 the result of filter 010

9、020030000.20.40.60.81 the signal of y 0123400.511.520123400.511.5(2). 三角序列（ 256 个点，两个周期）输入序列：低通滤波后的频谱和输出序列：0100200300-1-0.500.51 the result of filter 0100200300-1-0.500.51 the signal of y 0123405010015001234050100150(3). 矩形序列（两个周期，256个点，占空比0.5 ）输入序列：低通滤波后的频谱和输出序列：0100200300-2-1012 the result of filt

10、er 0100200300-1-0.500.51 the signal of y 0123405010015020001234050100150200主要代码：%y=ones(1,64) -ones(1,64) ones(1,64) -ones(1,64);%y=1 zeros(1,255);y=1:32 31:-1:-32 -31:32 31:-1:-32 -31:0./32; k=2*length(y); b,a=butter(3,0.2*pi); num1,den1=impinvar(b,a); h1,w=freqz(num1,den1); hh,tt=impz(b,a); yy1=co

11、nv(hh,y); %yy=filter(b,a,y);f=fft(y,k); ff1=fft(yy1,k); subplot(2,2,2); stem(yy1(1:length(y),.);title( the result of filter ); subplot(2,2,1); stem(y,.);title( the signal of y ); subplot(2,2,3); plot(w,abs(f); subplot(2,2,4); plot(w,abs(ff1); b.高通滤波：（1）矩形序列（ 256 个点，两个周期，占空比0.5 ）输入序列：高通滤波输出的频谱：高通滤波输出

12、：0100200300-1-0.500.51x 10-22 the result of filter 0100200300-1-0.500.51 the signal of y 012340501001502000123400.51x 10-20（2）矩形序列（256 个点，两个周期，占空比0.1 ）输入序列：高通滤波输出的频谱和高通滤波输出：0100200300-1-0.500.51x 10-22 the result of filter 0100200300-1-0.500.51 the signal of y 0123401002003000123400.511.5x 10-20主要代码

13、：y=ones(1,13) -ones(1,115) ones(1,13) -ones(1,115); k=2*length(y); b,a=butter(3,0.2*pi,high); num1,den1=impinvar(b,a); h1,w=freqz(num1,den1); hh,tt=impz(b,a,high); yy1=conv(hh,y); %yy=filter(b,a,y);f=fft(y,k); ff1=fft(yy1,k); subplot(2,2,2); stem(yy1(1:length(y),.);title( the result of filter ); sub

14、plot(2,2,1); stem(y,.);title( the signal of y ); subplot(2,2,3); plot(w,abs(f); subplot(2,2,4); plot(w,abs(ff1); 实验三一、实验名称fir 数字滤波器二、实验要求：设计一个截止频率为c=0.2的线性相位低通数字滤波器，1=0.3 ，2=0.3的线性相位带通滤波器，分别用矩形窗和海明窗对其进行截断，n为 61。输入序列64-128 点，输出 128-256 点。输入单位取样及矩形序列（占空比0.1 ） ，画出输出序列及其频谱。三、实验原理：理想低通数字滤波器，其频率特性为hd(ej)，

15、现假设其幅频特性|hd(ej)|=1，相频特性 ()=0，那么，该滤波器的单位抽样响应 hd(n)是以为 hd(0)为对称的 sinc 函数， hd(0)=c/。我们将 hd(n)截短，例如仅取hd(-m/2) ， hd(0)， hd(m/2) ，并将截短后的hd(n)移位，得h(n)=hd(n-m/2) n=0,1, ,m那么 h( n )是因果的，且为有限长，长度为m+1，令h(z)=h(n)z-n n=0,1,m 即得到所设计滤波器的转移函数。h(z) 的频率响应将近似hd(ej)，且是线形相位的。窗函数设计法是一种逼近，用其频响 h(ejw)去逼近所要求的理想滤波器频响hd(ejw),

16、用其有限长单位冲击响应h(n) 去逼近理想滤波器的无限长单位冲击响应hd(n) ，即：设计 fir df 的关键是求出h(n), 它应该是一个有限长因果序列。有限性可通过对hd(n) 截取一段， 即与某一窗函数相乘获得；因果性可通过在时域上进行的时延来获得，这不影响幅频特性，只影响相频。常用的窗函数有矩形窗、海明窗等。设计时，先根据算出 hd(n) ，再根据指定的窗函数点数以及窗的类型得出h(n)，对输入的待滤波序列和h(n)做卷积，即可达到滤波效果。具体实现时可根据线形卷积和圆周卷积的关系，通过补点把线形卷积化为圆周卷积，再根据离散时域的卷积定理，借助fft求出两序列的频谱，对其频域的乘积做

17、ifft, 即得到时域的圆周卷积。理想低通滤波器幅频特性可知：同理：带通滤波器的单位冲击响应为：h(n)=(sin 2n- sin 1n)/ ( n) 所以：10)()()(nmmnxmhny四、实验结果及分析：以下各图输入序列均为128 个点，输出序列均为256 个点，滤波器窗函数取样点的数目n 均为 61。低通滤波（1）单位冲击序列输入：02040608010012014000.51nx(n)输 入 序 列05010015020025030000.51n幅度/db加 矩 形 窗 低 通 滤 波 后 频 谱050100150200250300-0.200.2ny(n)加 矩 形 窗 低 通

18、滤 波 后 输 出 序 列02040608010012014000.51nx(n)输 入 序 列05010015020025030000.5n幅度/db加 海 明 窗 低 通 滤 波 后 频 谱050100150200250300-0.100.1ny(n)加 海 明 窗 低 通 滤 波 后 输 出 序 列00.10.20.30.40.50.60.70.80.91-30-20-100归 一 化 频 率 /p幅度/db加 矩 形 窗 后 响 应00.10.20.30.40.50.60.70.80.91-25-20-15-10-5归 一 化 频 率 /p幅度/db加 海 明 窗 后 响 应主要代码：

19、n=0:127;n=30;x(1)=1;x(2:128)=0; h=sin(0.2*pi*(n-30)./(pi*(n-30+eps);y=conv(x,h);y=abs(fft(y); subplot(3,1,1); stem(n,x); xlabel(n); ylabel(x(n); title(输入序列 ); subplot(3,1,2); stem(0:254,y) ;xlabel(n); ylabel( 幅度 /db); title( 加矩形窗低通滤波后频谱); subplot(3,1,3); stem(0:254,y); xlabel(n); ylabel(y(n); title(

20、加矩形窗低通滤波后输出序列); hd=sin(0.2*pi*(n-30)./(pi*(n-30+eps); w=0.5-0.5*cos(2*pi*n/128); h=hd.*w; y=conv(x,h); y=abs(fft(y); subplot(3,1,1) ;stem(n,x); xlabel(n); ylabel(x(n); title(输入序列 ); subplot(3,1,2); stem(0:254,y); xlabel(n); ylabel( 幅度 /db); title( 加海明窗低通滤波后频谱); subplot(3,1,3); stem(0:254,y); xlabel(

21、n); ylabel(y(n); title(加海明窗低通滤波后输出序列); n=0:127;h=sin(0.2*pi*(n-30)./(pi*(n-30+eps);h1,w1=freqz(h,1); subplot(2,1,1); plot(w1/pi,20*log10(abs(h1); xlabel(归一化频率 / ); ylabel( 幅度 /db); title( 加矩形窗后响应 ) d=sin(0.2*pi*(n-30)./(pi*(n-30+eps); w=0.5-0.5*cos(2*pi*n/128); h=d.*w; h1,w1=freqz(h,1); subplot(2,1,

22、2); plot(w1/pi,20*log10(abs(h1); xlabel(归一化频率 / ) ;ylabel( 幅度 /db); title( 加海明窗后响应 ) ; （2）矩形序列输入（占空比0.5）：020406080100120140-101nx(n)输 入 序 列050100150200250300050100n幅度/db加 矩 形 窗 低 通 滤 波 后 频 谱050100150200250300-202ny(n)加 矩 形 窗 低 通 滤 波 后 输 出 序 列020406080100120140-101nx(n)输 入 序 列05010015020025030002040n幅度/db加 海 明 窗 低 通 滤 波 后 频 谱050100150200250300-0.500.5ny(n)加 海 明 窗 低 通 滤 波 后 输 出 序 列主要代码：n=0:127;n=30;x=ones(1,32) -ones(1,32),ones(1,32) -ones(1,32); h=sin(0.2*pi*(n-n)./(pi*(n-n+eps);y=conv(x,h);y=abs(fft(y); subplot(3,1,1) ;stem(n,x); xlabel(n); ylabel(x(n); title( 输入序列 ); subplot(3,1,