求两个序列的卷积和。
编写仿真程序,并调试得到结果,进行分析。
程序:
①a=[1,0,1/3,0];%分母多项式系数
b=[1/6,1/2,1/2,1/6];%分子多项式系数
N=32;
n=0:
N-1;%一维数组,含N个分量,间隔为1
hn=impz(b,a,n);%单位冲激响应
gn=dstep(b,a,n);%单位阶跃响应
subplot(1,2,1),stem(n,hn,'k');
title('系统的单位冲激响应');
ylabel('h(n)');xlabel('n');
subplot(1,2,2),stem(n,gn,'k');
title('系统的单位阶跃响应');
ylabel('g(n)');xlabel('n');
②nf1=0:
20;
f1=0.8.^nf1;
subplot(2,2,1);stem(nf1,f1,'filled');
title('f1(n)');
nf2=0:
10;
lf2=length(nf2);
f2=ones(1,lf2);
subplot(2,2,2);stem(nf2,f2,'filled');
title('f2(n)');
y=conv(f1,f2);
subplot(2,1,2);stem(y,'filled');title('y(n)');
。
实验二频域采样
1.实验目的:
(1)掌握频率域采样会引起时域周期化的概念,以及频率域采样定理及其对频域采样点数选择的指导作用。
(2)会用MATLAB语言进行频域抽样与恢复时程序的编写方法。
2.实验原理:
了解频域采样定理的要点,掌握采样理论的结论:
“频域采样时域信号周期延拓”。
3.实验内容:
(1)频域采样理论的验证。
给定信号如下:
(2)编写程序分别对频谱函数
在区间
上等间隔采样32
和16点,得到
,再分别对
进行32点和16点IFFT,得到
。
(3)分别画出
、
的幅度谱,并绘图显示x(n)、
的波形,进行对比和分析,验证总结频域采样理论。
程序:
①
M=27;N=32;n=0:
M-1;
%产生M长三角波序列x(n)
xa=1:
ceil(M/2);%生成序列n+1
xb=floor(M/2):
-1:
1;%生成序列27-n;floor是向下取整,ceil是向上取整
xn=[xa,xb];
Xk=fft(xn,1024);%1024点FFT[x(n)],用于近似序列x(n)的FT
X32k=fft(xn,32);%32点FFT[x(n)]
x32n=ifft(X32k);%32点IFFT[X32(k)]得到x32(n)
X16k=X32k(1:
2:
N);%隔点抽取X32k得到X16(K)
x16n=ifft(X16k,N/2);%16点IFFT[X16(k)]得到x16(n)
subplot(3,2,2);
stem(n,xn,'.');%stem画离散序列图
boxon%给图形加边框
title('(b)三角波序列x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20])
k=0:
1023;wk=2*k/1024;%产生1024点DFT对应的采样点频率(关于π归一化值)
subplot(3,2,1);%画子图,3表示行数,2表示列数,1表示当前子图的序号数(以行元素优先顺序排列)。
plot(wk,abs(Xk));%绘制1024点DFT的幅频特性图
title('(a)FT[x(n)]');
xlabel('\omega/\pi');%'\'是转义符号等价于w/pi
ylabel('|X(e^j^\omega)|');%等价于e的jw次方
axis([0,1,0,200])%确定轴的范围,横轴从0~1,纵轴从0~到200
k=0:
N/2-1;
subplot(3,2,3);stem(k,abs(X16k),'.');boxon
title('(c)16点频域采样');xlabel('k');ylabel('|X_1_6(k)|');%下划线表示下标(转义字符‘_’)
axis([0,8,0,200])
n1=0:
N/2-1;
subplot(3,2,4);stem(n1,x16n,'.');boxon
title('(d)16点IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20])
k=0:
N-1;
subplot(3,2,5);stem(k,abs(X32k),'.');boxon
title('(e)32点频域采样');xlabel('k');ylabel('|X_3_2(k)|');axis([0,16,0,200])
n1=0:
N-1;
subplot(3,2,6);stem(n1,x32n,'.');boxon
title('(f)32点IDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])
②wp=0.25*pi;%滤波器的通带截止频率
ws=0.4*pi;%滤波器的阻带截止频率
Rp=1;As=15;%输入滤波器的通阻带衰减指标
ripple=10^(-Rp/20);Attn=10^(-As/20);
Fs=100;T=1/Fs;
Omgp=(2/T)*tan(wp/2);%原型通带频率预修正
Omgs=(2/T)*tan(ws/2);%原型阻带频率预修正
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s');%计算阶数n和3dB截止频率
[z0,p0,k0]=buttap(n);%归一化原型设计
ba=k0*real(poly(z0));%求原型滤波器系数b
aa=real(poly(p0));%求原型滤波器系数a
[ba1,aa1]=lp2lp(ba,aa,Omgc);%变换为模拟低通滤波器系数b,a
[bd,ad]=bilinear(ba1,aa1,Fs)%用双线性变换法求数字滤波器系数b,a
[sos,g]=tf2sos(bd,ad);%由直接型变换为级联型
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));%化为分贝值
subplot(2,2,1),plot(w/pi,abs(H));
ylabel('|H|');title('幅度响应');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/pi,angle(H)/pi);
ylabel('\phi');title('相位响应');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,3),plot(w/pi,dbH);title('幅度响应(dB)');
ylabel('dB');xlabel('频率(\pi)');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零极图');
4.思考题:
如果序列x(n)的长度为M,希望得到其频谱
在
上的N点等间隔采样,当N答:
先对原序列x(n)以N为周期进行周期延拓后取主值区序列,
再计算N点DFT则得到N点频域采样:
实验三信号的快速傅里叶变换实验
1、实验目的:
使学生进一步了解快速付里叶变换的理论;及用Matlab软件设计快速付里叶变换的方法,使学生进一步了解数字信号的频谱概念。
2、实验内容与步骤:
1)根据数字信号序列的长度决定2的整数次幂N,用Matlab软件设计出快速付里叶变换及反变换的程序,或按照范例程序进行修改,输入信号序列,运行程序,检查程序是否有错。
2)添加绘图语句,画出数字信号变换前后的波形。
运行程序,观察变换前后的波形有什么不同。
3)添加绘图语句,画出变换后的频谱图,运行程序,仔细观察频谱图,了解变换的意义。
3.实验程序及现象
m=10;
fori=1:
m+1//系统函数中B为分子多项式的系数,用矩阵的形式表示
ifi==1
B(i)=1;
elseifi==m+1
B(i)=-1;
elseB(i)=0;
end
end
A=[1,-1];//系统函数中A为分母多项式的系数
N=8192;
[H,f]=freqz(B,A,N);
plot(f*25/pi,abs(H));grid;
figure,plot(f*25/pi,angle(H));grid;
figure,zplane(B,A);
k=0:
N-1;f=2*k/N;
load('RawData.mat');
x=rawdata(1,1:
N);
w=filter(B,A,x);
X=abs(fft(x,N));
W=abs(fft(w,N));
figure;plot(x);title('输入信号');
figure;plot(w);title('输出信号');
figure;plot(f,abs(X));title('输入信号的幅频响应');
figure;plot(f,abs(W));title('输出信号的幅频响应');
figure;plot(f,abs(X),'b',f,abs(W),'r');
实验现象如下图
图形如下:
(3)极点分布图,一共有9个极点
(4)输入信号与(5)输出信号的不同,是由于输入信号经过了低通滤波器的滤波,滤除了不必要的波形。
(6)输入频响与(7)输出频响的区别,也与低通滤波器有关
快速傅里叶变换的目的:
快速傅里叶变化是为了为应用各种信号的实时处理提供条件的方法,大大的提高了傅里叶变化的计算量。
实验四用双线性变换法设计IIR数字滤波器
1.实验目的
(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。
(2)掌握用双线性变换法设计数字滤波器的计算机仿真方法。
(3)了解MATLAB有关双线性变换法的子函数。
2.实验涉及的MATLAB子函数
bilinear
功能:
双线性变换——将s域映射到z域的标准方法,将模拟滤波器变换成离散等效滤波器。
调用格式:
[numd,dend]=bilinear(num,den,Fs);将模拟域传递函数变换为数字域传递函数,Fs为取样频率。
3.实验原理
(1)双线性变换法是将整个s平面映射到整个z平面,其映射关系为
双线性变换法克服了脉冲响应不变法从s平面到z平面的多值映射的缺点,消除了频谱混叠现象。
但其在变换过程中产生了非线性的畸变,在设计IIR数字滤波器的过程中需要进行一定的预修正。
(2)双线性变换法设计IIR数字滤波器的步骤:
①输入给定的数字滤波器设计指标;
②根据公式进行预修正,将数字滤波器指标转换成模拟滤波器设计指标;
③确定模拟滤波器的最小阶数和截止频率;
④计算模拟低通原型滤波器的系统传递函数;
⑤利用模拟域频率变换法,求解实际模拟滤波器的系统传递函数;
⑥用双线性变换法将模拟滤波器转换为数字滤波器。
4.实验内容
采用双线性变换法设计一个巴特沃斯数字低通滤波器,要求:
滤波器采样频率
5.思考题
用双线性变换法设计数字滤波器过程中,变换公式中T的取值,对设计结果有无影响?
为什么?
答:
没有,一般取2/T=1,方便计算
6.实验程序现象
程序如下:
wp=0.25*pi;%滤波器的通带截止频率
ws=0.4*pi;%滤波器的阻带截止频率
Rp=1;As=15;%输入滤波器的通阻带衰减指标
ripple=10^(-Rp/20);Attn=10^(-As/20);
Fs=100;T=1/Fs;
Omgp=(2/T)*tan(wp/2);%原型通带频率预修正
Omgs=(2/T)*tan(ws/2);%原型阻带频率预修正
[n,Omgc]=butter(Omgp,Omgs,Rp,As,'s');%计算阶数n和3dB截止频率
[z0,p0,k0]=buttap(n);%归一化原型设计
ba=k0*real(poly(z0));%求原型滤波器系数b
aa=real(poly(p0));%求原型滤波器系数a
[ba1,aa1]=lp2lp(ba,aa,Omgc);%变换为模拟低通滤波器系数b,a
[bd,ad]=bilinear(ba1,aa1,Fs)%用双线性变换法求数字滤波器系数b,a
[sos,g]=tf2sos(bd,ad);%由直接型变换为级联型
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));%化为分贝值
subplot(2,2,1),plot(w/pi,abs(H));
ylabel('|H|');title('幅度响应');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/pi,angle(H)/pi);
ylabel('\phi');title('相位响应');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,3),plot(w/pi,dbH);title('幅度响应(dB)');
ylabel('dB');xlabel('频率(\pi)');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零极图');
由频率特性曲线可知,该设计结果在通阻带截止频率处能满足Rp≤1 dB、As≥15 dB的设计指标要求,系统的极点全部在单位圆内,是一个稳定的系统。
由n=5可知,设计的巴特沃斯数字低通滤波器是一个5阶的系统,原型Ha(s)在s=-∞处有5个零点,映射到z=-1处。
双线性变换法特点:
双线性变换法克服了脉冲响应不变法从s平面到z平面的多值映射的缺点,消除了频谱混叠现象。
但其在变换过程中产生了非线性的畸变,在设计IIR数字滤波器的过程中需要进行一定的预修正。
总结
第一次做DSP的实验,对matlab软件非常陌生,其程序也是一窍不通,但是通过老师对软件使用方法的讲解和源代码的开放,使我们对数字信号处理实验有了初步的认识,开始的时候始终不能理解程序,用了老师告诉我们的设置断点单步运行的方法后,才对程序有了一定的了解。
通过本次实验熟悉掌握了求系统响应的方法,在时域中,已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应。
系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。
还有频域采样原理,频域采样点数N必须大于等于时域离散信号的长度M(即N≥M),才能使时域不产生混叠。
并且可知“频域采样时域信号周期延拓”。
第二次实验,老师要我们补充一个不完整的程序,实现信号的快速傅里叶变换,开始的时候总是不明白我们需要补充的是什么内容,后来通过老师的一番讲解,才知道通过一段循环程序写出一个系统函数的表达式,之前以为这个式子根本可以不用循环就能实现,直接用矩阵写出,后来明白这种方法不具有普遍性,一旦遇到阶数比较多的函数这个方法就行不通了。
通过两次实验,熟悉了matlab7.0的使用方法,以及函数的建立和绘图。
进一步了解了DFT和FFT,以及二者的误差。
并掌握了如何从一个信号中获得清楚的谱线。
在整个实验过程中,书本上的知识得到进一步巩固,熟悉了matlab的画图技巧,为以后的学习和实验打下了基础。