《测试信号分析与处理》
实验一 差分方程、卷积、z变换
一、 实验目的
通过该实验熟悉 matlab软件的基本操作指令,掌握 matlab软件的使用方法,掌握数字信号处理中的基本原理、方法以及matlab函数的调用。
二、 实验设备
1、微型计算机1台; 2、matlab软件1套
三、 实验原理
Matlab 软件是由mathworks公司于1984年推出的一套科学计算软件,分为
总包和若干个工具箱,其中包含用于信号分析与处理的sptool工具箱和用于滤波器设计的fdatool工具箱。它具有强大的矩阵计算和数据可视化能力,是广泛应用于 信号分析与处理中的功能强大且使用简单方便的成熟软件。Matlab软件中已有大量的关于数字信号处理的运算函数可供调用,本实验主要是针对数字信号处理中的差分方程、卷积、z变换等基本运算的matlab函数的熟悉和应用。 差分方程(difference equation)可用来描述线性时不变、因果数字滤波器。用x表示滤波器的输入,用y表示滤波器的输出。
a0y[n]+a1y[n-1]+…+aNy[n-N]=b0x[n]+b1x[n-1]+…+bMx[n-M] (1)
ak,bk 为权系数,称为滤波器系数。
N为所需过去输出的个数,M 为所需输入的个数
卷积是滤波器另一种实现方法。
y[n]= ∑ x[k] h[n-k] = x[n]*h[n] (2)
等式定义了数字卷积,*是卷积运算符。输出 y[n] 取决于输入 x[n] 和系统的脉冲响应h[n]。
传输函数H(z)是滤波器的第三种实现方法。
H(z)=输出/输入= Y(z)/X(z) (3) 即分别对滤波器的输入和输出信号求z变换 ,二者的比值就是数字滤波器的传输函数。
序列x[n]的z变换定义为
X (z)=∑x[n]z-n (4)
把序列 x[n] 的 z 变换记为Z{x[n]} = X(z)。 由 X(z) 计算 x[n] 进行 z 的逆变换 x[n] = Z-1{X(z)}。
Z 变换是 Z-1 的幂级数,只有当此级数收敛,Z 变换才有意义,而且同
一个 Z 变换等式,收敛域不同,可以代表不同序列的 Z 变换函数。
这三种数字滤波器的表示方法之间可以进行相互转换。
四、 实验步骤
1、熟悉matlab软件基本操作指令。读懂下列matlab程序指令,键入程序并运行,观察运行结果。
Conv.m% 计算两个序列的线性卷积;
%----------------------------------------------------------------- clear;
N=5; M=6;
L=N+M-1; x=[1,2,3,4,5]; h=[6,2,3,6,4,2]; y=conv(x,h); nx=0:N-1; nh=0:M-1; ny=0:L-1;
subplot(231);
stem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on; subplot(232);
stem(nh,h,'.k');xlabel('n');ylabel('h(n)');grid on; subplot(233);
stem(ny,y,'.k');xlabel('n');ylabel('y(n)');grid on;
filter.m;%求一个离散系统的输出; clear;
x=ones(100); t=1:100;
b=[.001836,.007344,.011016,.007374,.001836]; a=[1,-3.0544,3.8291,-2.2925,.55075]; y=filter(b,a,x); clear;
impz .m% 计算滤波器的冲击响应
b=[.001836,.007344,.011016,.007374,.001836]; a=[1,-3.0544,3.8291,-2.2925,.55075]; [h,t]=impz(b,a,40); subplot(221)
stem(t,h,'.');grid on; ylabel('h(n)')
xlabel('n')
filter.m% 计算滤波器的阶跃响应 x=ones(100);t=1:100; y=filter(b,a,x); subplot(222)
plot(t,x,'g.',t,y,'k-');grid on; ylabel('x(n) and y(n)') xlabel('n')
例题运行结果图
66806044h(n)x(n)y(n)200n5402002002n405n10
2、编程求出下列问题的解 1)、滤波器的差分方程为:y[n]=x[n]-0.8x[n-1]-0.5y[n-1]
求出此滤波器脉冲响应和阶跃响应的前十个采样值。
clear;
%impz.m% 计算滤波器的冲击响应 b=[1,-.8]; a=[1,.5];
[h,t]=impz(b,a,10); stem(t,h,'.');gird on; ylabel('h(n)')
xlabel('n')
10.50h(n)-0.5-1-1.501234n567clear;
%filter.m% 计算滤波器的阶跃响应 x=ones(10);t=1:10; b=[1,-.8]; a=[1,.5];
y=filter(b,a,x);
plot(t,x,'g.',t,y,'k.');gird on; ylabel('x(n) and y(n)') xlabel('n')
10.80.6x(n) and y(n)0.40.20-0.2-0.4123-n
45n6710 2)、系统的脉冲响应为h[n]=e(u[n]-u[n-3]),用卷积求系统的阶跃响应。 N=25; M=3;
L=N+M-1;
x=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]; h=[1,.3679,.1353]; y=conv(x,h); nx=0:N-1; nh=0:M-1; ny=0:L-1; subplot(231);
stem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on; subplot(232);
stem(nh,h,'.k');xlabel('n');ylabel('h(n)');grid on; subplot(233);
stem(ny,y,'.k');xlabel('n');ylabel('y(n)');grid on;
1121.5h(n)x(n)0.50.5y(n)10.50020n40001n20020n4
五、实验讨论和分析
1、差分方程、卷积、z变换和傅里叶变换之间如何进行转换?
答:差分方程;a0y[n]+a1y[n-1]+a2y[n-2]+`````+aNy[n-N]=b0x[n]+b1x[n-1]+……+bMx[n-M] 卷积是由输入x[n]所引起的全部输出y[n]是所有这些加权脉冲相应之和。即y{n}=x[n]*h[n]只要知道脉冲响应和输入就可以得到输出 Z变换是把时域信号向频域进行转换X(z)=∑x[n]zˇ-n Y(z)=∑y[n]zˇ-n 脉冲响应是传输函数的逆z变换 傅里叶变换X(Ω)=∑x[n]eˇ-jnΩ 2、边界效应是如何产生的?它对信号的滤波效果有何影响?
答:多数情况下,采样开始之前的输入情况是未知的,当脉冲响应与未知的的输
入采样点重叠时,由于实际的输出值可能受采样开始之前输入信号的影响,所以无法准确的计算输出。计算的开始和末尾都存在这种现象。仅当输入序列与脉冲响应完全重叠时,计算才有意义,这种现象就是边界效应。当一个系统开始运行或条件改变时,输出需要一些时间过渡到新的稳态。 边界效应会产生输出的暂态部分和稳态部分,会影响滤波效果,并且会导致失真现象出现。
实验二 数字滤波器综合设计
一、 实验目的
通过该设计实验掌数字滤波器设计的一般步骤,掌握利用matlab 软件设计数字滤波器的方法,熟悉sptool工具箱的使用方法。
二、 实验设备
1、微型计算机1台; 2、matlab软件1套
三、 实验原理
一)、滤波器的形状及重要参数
理想滤波器的形状是矩形,图 1 给出非理想滤波器。
图 1
通带:增益高的频率范围,信号可以通过,称为滤波器的通带。
阻带:增益低的频率范围,滤波器对信号有衰减或阻塞作用,称滤波器的阻带。 滤波器截止频率:增益为最大值的0.707倍时所对应的频率为滤波器截止频率 增益通常用分贝(dB)表示。增益(dB)= 20log(增益) 增益为 0.707 时对应 -3dB,因此截止频率常被称为 -3dB。 滤波器的带宽:对于低通滤波器宽带是从0 ~ - 3dB
对于高通滤波器宽带是从 - 3dB~采样频率的一半
对于带通滤波器带宽是截止频率之间的频率距离
二)加窗低通 FIR 滤波器的设计
1. 在过渡带宽度的中间,选择通带边缘频率(Hz): f1=所要求的通带边缘频率+(过渡带宽度)/2
2. 计算 Ω1=2πf1/fs,并将此值代入理想低通滤波器的脉冲响应 h1[n] 中: h1[n] = sin(nΩ1)/nπ
3. 从表中选择满足阻带衰减及其他滤波器要求的窗函数,用表中 N 的公式计算所需要的非零项数目。选择奇数项,这样脉冲响应可以完全对称,避免了滤波器产生相位失真,对于|n|≤(N-1)/2,计算窗函数w[n]。
4. 对于|n|≤(N-1)/2,从式 h[n]=h1[n]w[n]计算(有限)脉冲响应,对于其他 n 值h[n]=0,此脉冲响应是非因果的。
5. 将脉冲响应右移 (N-1)/2,确保第一个非零值在n=0处,使此低通滤波器为因果的。 三)、设计低通巴特沃斯滤波器:
1) 确定待求通带边缘频率 fp1 Hz 、待求阻带边缘频率fs1 Hz 和待求阻带衰减 - 20logδsdB(或待求阻带增益 20logδsdB)。通带边缘频率对应 –3dB增益。
2) 用式 Ω=2πf/fs 把由 Hz 表示的待求边缘频率转成由弧度表示的数字频率,得到 Ωp1 和Ωs1 。
3) 计算预扭曲模拟频率以避免双线性变化带来的失真。由 ω=2fs tan(Ω/2) 求得 ωp1和 ωs1,单位是弧度/秒。
4) 由已给定的阻带衰减 - 20logδs(或增益- 20logδs)确定阻带边缘增益 δs 。
5) 计算所需滤波器的阶数n 取整数。
6)把 ωp1代入 n 阶模拟巴特沃斯滤波器传输函数H(s)中,并对 H(s) 进行双线性变换得到 n 阶数字传输函数 H(z)。滤波器实现所需的差分方程可直接从传输函数 H(s) 求出。 。 四)、低通切比雪夫Ⅰ型滤波器的设计:
1)确定待求的通带与阻带边缘频率 fp1 和fs1 、待求的通带边缘增益 20log(1- δp) 和待求的阻带衰减-20logδs(或待求的阻带增益 20logδs )。 2)用公式 Ω=2πf/fs 将待求的边缘频率转换为数字频率(用弧度表示),得到 Ωp1 和 Ωs1 。 3)对数字频率采用预扭曲以避免双线性变换引起的误差。由 ω=2fs tan(Ω/2) 得到ωp1和 ωs1,单位是弧度/秒。
4)由指定的通带边缘增益 20log(1- δp) ,确定通带边缘增益 1- δp 。计
算参数ε。
5)由指定的衰减-20logδs(或增益 20logδs),确定阻带边缘增益 δs 。 6)计算所需的阶数n。
7)将 ωp1 和 δp 代入 n 阶模拟切比雪夫Ⅰ型滤波器的传输函数 H(s),并对其进行双线性变换,得到 n 阶数字滤波器传输函数 H(z)。实现滤波器所需的差分方程可由传输函数 H(z) 直接得到。
四、实验步骤
1、
任选第9、10章后滤波器设计题各2题,利用matlab编程完成滤波器的设计,并画出滤波器的脉冲响应、幅度响应和相位响应图。
习题9.15
f1=4000;%信号频率Hz f2=5000;%信号频率Hz f3=6000;%信号频率Hz fs=12000;%采样频率Hz N=32;%采样点数
t=(0:N-1)/fs;%采样时间
x1=sin(2*pi*f1*t);%信号采样值 x2=sin(2*pi*f2*t);%信号采样值 x3=sin(2*pi*f3*t);%信号采样值 x=x1+x2+x3; y=filter(h,1,x);
f1=3000+250; fs=12000; w=2*f1/fs; n=3.32*fs/500; h=makelp(n,w,'hanning'); [mag,phase,w]=dtft(h); plot(t,x,'g',t,y,'k-')
老师,这道题的错误不会改,不能运行。
习题9.23:
h=bandfilt(59,0.31875,0.68125,1,'hanning'); [mag,phase,w]=dtft(h); plotdtft(mag,phase,w,2); stem(0:116,h,'.'); ylabel('h(n)'); xlabel('n');
习题10.6
n=buttord(0.25,0.375,3,44); [b,a]=butter(n,0.25);
[mag,phase,w] = dtft(b,a); plotdtft(mag,phase,w,1); [h,t]=impz(b,a,40); subplot(111)
stem(t,h,'.');grid on; ylabel('h(n)') xlabel('n')
习题10.12
n=buttord(0.25,0.375,3,44); [b,a]=butter(n,0.25);
[mag,phase,w] = dtft(b,a); plotdtft(mag,phase,w,1); [h,t]=impz(b,a,40); subplot(111)
stem(t,h,'.');grid on; ylabel('h(n)') xlabel('n')
五、实验讨论和分析
1、设计得到的滤波器与设计要求有无差别?如果有,请分析误差产生的原因。 答:有差别。在设计FIR滤波器时,我们不可能得到理想的滤波器,而是要选用合适的窗函数,来满足阻带衰减要求,加窗后滤波器形状就不是理想的了,并且在它的通带和阻带内有波纹,还有就是滤波器系数自身的量化,如果选用比特数少,就会产生大的误差,量化也会影响IIR的稳定性,IIR滤波器不能保证无相位失真
2、 FIR滤波器与IIR滤波器的优缺点分别是什么?针对具体信号进行滤波时,如何选择?
答:FIR滤波器的最主要的特点是没有反馈回路,故不存在不稳定的问题;同时,可以在幅度特性是随意设置的同时,保证精确的线性相位。稳定和线性相位特性是FIR滤波器的突出优点。另外,它还有以下特点:设计方式是线性的;硬件容易实现;滤波器过渡过程具有有限区间;相对IIR滤波器而言,阶次较高,其延迟也要比同样性能IIR滤波器大得多。IIR滤波器的首要优点是可在相同阶数时取得更好的滤波效果。但是IIR滤波器设计方法的一个缺点是无法控制滤波器的相位特性。由于极点会杂散到稳定区域之外,自适应IIR滤波器设计中碰到的一个大问题是滤波器可能不稳定。因此,一般采用FIR滤波器作为自适应滤波器的结构。
实验三 数字信号处理综合设计 一.对实际信号处理
1.语音信号的频谱分析
要求首先画出语音信号的时域波形;然后对语音信号进行频谱分析,在MATLAB中,
可以利用函数fft对信号进行快速付立叶变换,得到信号的频谱特性;从而加深对频谱特性的理解。
fs=20000; %语音信号采样频率为20000
x1=wavread('d:\\lianxi.wav',20000); %读取语音信号的数据,赋给变量x1 sound(x1,20000); %播放语音信号
y1=fft(x1,1024); %对信号做1024点FFT变换 f=fs*(0:511)/1024; figure(1)
plot(x1) %做原始语音信号的时域图形 title('原始语音信号'); xlabel('time n'); ylabel('fuzhi n'); figure(2)
freqz(x1) %绘制原始语音信号的频率响应图 title('频率响应图') figure(3) subplot(2,1,1);
plot(abs(y1(1:512))) %做原始语音信号的FFT频谱图 title('原始语音信号FFT频谱') subplot(2,1,2); plot(f,abs(y1(1:512))); title('原始语音信号频谱') xlabel('Hz'); ylabel('fuzhi');
2. 以低通滤波器为例,对信号进行处理,回放语音信号在MATLAB中,函数sound可以对声音进行回放。其调用格式: sound(x,fs,bits);
可以感觉滤波前后的声音有变化。
低通:
fs=20000;
x1=wavread('d:\\lianxi.wav',20000); t=0:1/20000:(size(x1)-1)/20000; wp=0.1*pi; ws=0.5673*pi; Rp=1; Rs=100; Fs=20000; Ts=1/Fs;
wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标 ws1=2/Ts*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s'); %选择滤波器的最小阶数 [Z,P,K]=buttap(N); %创建butterworth模拟滤波器 [Bap,Aap]=zp2tf(Z,P,K); [b,a]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 [H,W]=freqz(bz,az); %绘制频率响应曲线 figure(1)
plot(W*Fs/(2*pi),abs(H)) grid
xlabel('频率/Hz') ylabel('频率响应幅度') title('Butterworth') f1=filter(bz,az,x1); figure(2)
subplot(2,1,1)
plot(t,x1) %画出滤波前的时域图 title('滤波前的时域波形'); subplot(2,1,2)
plot(t,f1); %画出滤波后的时域图 title('滤波后的时域波形');
sound(f1,20000); F0=fft(f1,1024); f=fs*(0:511)/1024; figure(3) y1=fft(x1,1024); subplot(2,1,1);
plot(f,abs(y1(1:512))); title('滤波前的频谱') xlabel('Hz'); ylabel('fuzhi'); subplot(2,1,2)
F1=plot(f,abs(F0(1:512))); title('滤波后的频谱') xlabel('Hz'); ylabel('fuzhi');
%播放滤波后的信号 %画出滤波前的频谱图 %画出滤波后的频谱图
声音略
2分析结果
答:这次综合设计,熟悉了Matlab的用法,对于数字信号处理有了更加深刻的认识。对于不同的滤波器,FIR滤波器,如窗函数滤波器;IIR滤波器,如切比雪夫I型滤波器,巴特沃兹滤波器都有了更加形象的认识。在设计的过程中也遇到了很多的问题,比如不知道应该采用什么样的函数,函数的用法,以及如何根据要求画出合格的波形。为了更加形象准确,滤波器波形的横坐标应该是频率值,而纵坐标应该是增益(dB),才能使滤波器的特性曲线能够充分反映出滤波器的特点。
因篇幅问题不能全部显示,请点此查看更多更全内容