陆宁-201105090108-数字信号处理实验报告

数字信号处理实验报告姓名：陆宁学号：201105090108实验一信号、系统及系统响应题目：当系统的输入输出差分方程为：通过MATLAB编程实现并考虑如下问题：（1）当系统的输入为单位冲激函数时，分别利用filter函数和impz函数得到的系统单位冲激响应曲线。（2）当系统的输入为单位阶跃函数时，分别利用filter函数和impz函数得到的系统单位阶跃响应曲线。（3）对于不同的输入，系统的输出有什么变化，试讨论之。要求：提交报告，报告包括：原理、程序、结果及分析。原理：系统的零状态响应就是在系统初始状态为零的条件下微分方程的解。在MATLAB中，可以用函数y=filter(p,d,x)实现差分方程的仿真，其中p,d分别为系统差分方程左右端系数向量，x为数输入向量；impz函数可直接求出离散系统的单位冲激响应，为：impz(p,d)。程序：（1）用filter函数a1=[1,-0.8,-0.5];b1=[0.7,0.3];n=0:20;x1=[1zeros(1,20)];y1filter=filter(b1,a1,x1);stem(n,y1filter);title('y1filter');xlabel('x');ylabel('y');用impz函数a1=[1,-0.8,-0.5];b1=[0.7,0.3];impz(b1,a1,21);(2)用filter函数a1=[1,-0.8,-0.5];b1=[0.7,0.3];x2=ones(1,21);y1filter=filter(b1,a1,x2);stem(n,y1filter);title('y1filter_step');xlabel('x');ylabel('y');用impz函数a1=[1,-0.8,-0.5];b1=1;impz(b1,a1);（3）即y=filter(p,d,x)用来实现差分方程，d表示差分方程输出y的系数，p表示输入x的系数，而x表示输入序列。输出结果长度数等于x的长度。而y=impz(p,d,N)是用来实现冲击响应的，d和p的定义见filter，N表示冲击响应输出的序列个数。结果：（1）用filter函数用impz函数（2）用filter函数用impz函数分析：输入为离散信号时，输出为连续；相反，输入连续信号时，输出为离散。对于单位阶跃和单位抽样而言，输出没有变化。实验二使用FFT作谱分析题目：编程实现下列问题并讨论所得到的结果。（1）使用FFT对MATLAB中的数据noissin.dat进行谱分析。（2）使用FFT对语音数据noisyspeech.wav进行谱分析。要求：提交报告，报告包括：原理、程序、结果及分析。原理：（1）离散傅里叶变换(DFT)公式为：X(k)=∑x(n)W^nk;x(n)=∑X(k)W^-nk;其中w=e^(2∏nk/N),N为离散序列的长度。（2）快速傅里叶变换(FFT)是利用w因子的取值特点，减少DFT的复数乘法的次数。其中一种是时间抽取基2算法，它将时间按奇偶逐级分开，直到两点的DFT。MATLAB提供了fft函数可用于计算FFT，器调用形式为;X=fft（x）或X=fft（x，N），N为2的整数次幂，若x的长度小于N，则补零，若超过则舍去N后的数据。（3）函数形式[y,fs,bits]=wavread('Blip',[N1N2]);用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。[N1N2]表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值）。函数形式s=noissin（n1，n2）用于读取MALAB的噪声信号。注：noissin.dat含有噪声的正弦信号程序：(1)用FFT对noissin.data做频谱分析loadnoissin;s=noissin(1:50);S=fft(s);subplot(2,1,1);plot(abs(s));xlabel('n');ylabel('|s(n)|');gridon;subplot(2,1,2);plot(abs(S));xlabel('k');ylabel('|S(k)|');gridon;(2)用FFT对noisyspeech.wav进行频谱分析xx=wavread('noisyspeech.wav');fs=100;N=128;x=xx(1:N);n=1:N;X=abs(fft(x,N));subplot(2,1,1);plot(n,x);xlabel('n');ylabel('|x(n)|');gridon;subplot(2,1,2);plot(n,X);xlabel('k');ylabel('|X(k)|');gridon;结果：(1)用FFT对noissin.data做频谱分析(2)用FFT对noisyspeech.wav进行频谱分析分析：经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是2л／N≤D。可以根据此式选择FFT的变换区间N。误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时，离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。对模拟信号进行频谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的普分析进行。实验三使用双线性Z变换设计IIR滤波器题目：使用双线性Z变换法设计一个低通数字IIR滤波器，给定的数字滤波器的技术指标为，抽样频率为。要求：提交报告，报告包括：原理、程序、结果及分析。原理：（1）设计滤波器就是要设计一个系统是其能让一定频率的波段通过或滤去，对IIR滤波器，器转移函数是：H(Z)=(∑bz^(-r))/(1+∑az^(-k))。(2）设计的一般原则：若使滤波器拒绝某个频率，应在单位园上相应的频率处设置一个零点，反之则设置一个极点。(3)低通数字IIR滤波器设计步骤：给出数字低通滤波器的技术要求→映射为模拟低通的技术要求→归一化为模拟低通滤波器的技术要求→设计出G(P)→G(S)→映射到数字滤波器的转移函数G(Z)。(4)双线性Z变换，即S平面到Z平面的映射关系：S=2（Z-1）/Ts（Z+1）。程序：fp=100;fs=300;Fs=1000;rp=3;rs=20;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs=Fs/Fs;wap=tan(wp/2);was=tan(ws/2);[n,wn]=buttord(wap,was,rp,rs,'s');[z,p,k]=buttap(n);[bp,ap]=zp2tf(z,p,k);[bs,as]=lp2lp(bp,ap,wap);[bz,az]=bilinear(bs,as,Fs/2);[h,w]=freqz(bz,az,256,Fs*1000);plot(w,abs(h));xlabel('w');ylabel('|h(w)|');gridon;bpapbsasbzaz结果：bp=001ap=1.00001.41421.0000bs=0.1056as=1.00000.45950.1056bz=0.06750.13490.0675az=1.0000-1.14300.4128分析：(1)[n,wn]=buttord(wap,was,rp,rs,'s');%用于计算巴特奥斯数字滤波器的阶数n和3dB截止频率wc。调用参数wap，was分别为数字滤波器的通带、阻带截止频率的归一化值，要求：0≤wap≤1，0≤was≤1。当was≤wap时，为高通滤波器；n，wn作为butter函数的调用参数。(2)bilinear为双线性插值。在IIR数字滤波器设计中有把模拟的原型滤波器转变成相应的数字滤波器，即通过已知S变换中模拟滤波器的系数求数字滤波器的系数。其中有两种方法来进行这种变换，一种是激冲不变法，另一种是双线性Z变换法。bilinear函数便是双线性Z变换法。函数中的NUMs和DENs是原模拟滤波器S变换中的分子和分母系数，A和B是数字滤波器的系数。实验四使用窗函数发设计FIR滤波器题目：根据下列指标，设计一个FIR数字低通滤波器：ω=2.0∏，ω=4.0∏Ap=0.25dB，As=50dB。（1）分别考查矩形窗、三角窗、Hanning窗、Hamming窗并分析这些不同的窗函数对滤波器的幅度响应有什么影响。（2）选择一个合适的窗函数，确定单位冲激响应，绘出所设计的滤波器的幅度响应。原理：（1）窗函数法：将理想低通数字滤波器Hd（z）的单位抽样响应hd（n）截短并位移，可得到对理想滤波器的近似的转移函数H（z）。此过程即为在hd(n)上施加矩形窗口。（2）FIRDF设计的窗含函数法仅能给出通带截止频率wp，其他几个参数是靠h（n）的长度及所使用的窗含数的性能来决定的。常用的窗函数有：矩形窗、三角窗、汉宁窗、汉明窗。MATLAB窗函数法设计FIR数字低通滤波器用函数firl其调用形式为：b=firl（N，Wn）。式中N为滤波器的阶次，Wn是通带截止频率，b是设计好的滤波器系数h（n）。（3）freqz函数，用来求系统的频率响应，其调用形式为：[H,w]=freqz(b,a,N,’whole’,Fs)。常用的窗函数：1、矩形窗函数为boxcar和rectwin，调用格式：w=boxcar（N）w=rectwin(N)其中N是窗函数的长度，返回值w是一个N阶的向量。2、三角窗函数为triang，调用格式：w=triang(N)3、汉宁窗函数为hann，调用格式：w=hann（N）4、海明窗函数为hamming，调用格式：w=hamming（N）程序：（1）生成四种窗函数矩形窗、三角窗、汉宁窗、海明窗，并观察其频率响应。（1）矩形窗n=20;window=rectwin(n);[h,w]=freqz(window,1);subplot(4,2,1);stem(window);title('矩形窗');subplot(4,2,2);plot(w/pi,20*log(abs(h))/abs(h(1)));title('矩形窗的频率响应');holdon;（2）三角窗window=triang(n);[h,w]=freqz(window,1);subplot(4,2,3);stem(window);title('三角窗');subplot(4,2,4);plot(w/pi,20*log(abs(h))/abs(h(1)));title('三角窗的频率响应');holdon;（3）hanningwindow=hann(n);[h,w]=freqz(window,1);subplot(4,2,5);stem(window);title('汉宁窗');subplot(4,2,6);plot(w/pi,20*log(abs(h))/abs(h(1)));title('汉宁窗的频率响应');holdon;(4)hammingwindow=hamming(n);[h,w]=freqz(window,1);subplot(4,2,7);stem(window);title('海明窗');subplot(4,2,8);plot(w/pi,20*log(abs(h))/abs(h(1)));title('海明窗的频率响应');holdon;结果：分析：通过对窗函数的调用，使4个窗函数的时域及频域图放置在一起，比较直观地看出各个窗函数的特性。（2）根据下列技术指标，设计一个FIR数字低通滤波器：wp=0.2π，ws=0.4π，ap=0.25dB，as=50dB，选择一个适当的窗函数，确定单位冲激响应，绘出所设计的滤波器的幅度响应。提示：根据窗函数最