数字信号处理 MATLAB 上机实验 教学讲义

74
4.3 离散卷积
• matlab 中卷积运算的指令是c=conv(a,b),其中a和b是有 限长的序列。
【例】计算下面两个序列的卷积
a=[-2 0 1 –1 3]; b=[1 2 0 -1]; 两个序列的都从0开始的
75
%conva.m x=[-2 0 1 -1 3]; h=[1 2 0 -1 0]; y=conv(x,h); N=length(y)-1; n=0:1:N; stem(n,y,'fill'); grid on; title('conv of x and h'); ylabel('y'); xlabel('Time index n');
/6
的正
71
• 单位脉冲序列 (n) 和单位阶跃序列u(n) 可以用ones(1,n)和zeros(1,n)来生成单位脉冲 序列和单位阶跃序列 ones(1,n)命令产生1行n列的1值 zeros(1,n)命令产生1行n列的0值
72
【例】 产生单位脉冲序列 (n 3) 和单位阶跃序列u(n-3)
1
1.5 2 frequency in rad
2.5
3
95
【例】比较高密度谱和高分辨率谱之间的差异，考虑
x（n）=cos(0.48*pi*n)+cos(0.52*pi*n)
（１）取x（n）（0 n 10 ）时，求x（n）的DFT。 （２）将（１）中的x（n）以补零方式使x（n） 加 长到 0 n 100 ，求x（n）的DFT 。 （３）取x（n）（ 0 n 100 ）时，求x（n）的 DFT画出幅频特性。
Matlab基础介绍
一 二 三 四 、Matlab 简介 、Matlab 的安装与启动 、Matlab 编程基础 、Matlab 在数字信号处理课 程中的应用
1
四、 Matlab 在数字信号处理课程中的应用 举例 • • • • 常见连续信号 离散信号 离散卷积 LTI系统
65
4.1 常见连续信号
n=-2:30; x=[zeros(1,5),1,zeros(1,27)]; y=[zeros(1,5),ones(1,28)]; subplot(2,1,1); stem(n,x,'fill');grid on; subplot(2,1,2) stem(n,y,'fill');grid on;
73
68
69
• 随机信号：
【例】 rand产生均匀分布的白嘈声，randn产生高斯分
布的白噪声 t=0:0.01:1; y=randn(1,length(t)); plot(t,y);grid on;
70
4.2离散信号
•
周期序列 弦信号
【例】产生一个幅度为2，频率为4Hz,相位为 A=2; phi=pi/6; omega=2*pi/12; n=-10:10; x=A*sin(omega*n+phi); stem(n,x,'fill'); grid on;
stem(t,h2,'fill');
grid on; title('单位冲激响应2');
ylabel('h2(n)');
xlabel('n')
87
figure(2) zplane(num,den); %作出极－零图
figure(3) freqz(num,den,N,Fs); grid on; title('频率响应'); %做出幅频和相频响应曲线
76
77
【例】有序列
1 A(n) 0
1 B(n) 0
n 3,4,,12 else
n 2,3,,9 else
和
百度文库
求这两个完整序列的卷积，并图示。
78
%完整序列卷积%conv1.m a=ones(1,10); n1=3;n2=12; %完整a(n)序列的非平凡值和区间端点
88
89
90
91
【例】有一调幅信号
xa (t ) [1 cos(2 100t )]cos(2 600t )
•用做频谱分析，要求能分辨xa(t)的所有频率分量， •若用fs=3kHz频率抽样，抽样数据为５１２点， •做频谱分析，求X(k)=DFT[x(n)]，５１２点， •画出X(k)的幅频特性|X(k)|。
82
freqz：
因果LTI系统的的幅频响应和相频响应曲线，在matlab 中可以用以下函数作出： [H,w] =freqz(num,den,N)，频率响应在单位圆的上半周取N 个等分点； [H,w]=freqz(num,den, N,'whole')，频率响应在单位圆取N 个等分点； H=freqz(num,den, w) ，频率响应在w设定的频率上。 该函数可以同时作出幅频和相频响应图 其中： num-----系统传递函数分子系数组成的行向量； den------系统传递函数分母系数组成的行向量； N ---------是频率响应的点数，最好为2的幂缺省值512； w --------是返回频率轴坐标向量。
97
% part 2
% High density spectrum (100 samples) based on the first 10 samples of x(n) figure(2)
n3=[0:1:99];x=cos(0.48*pi*n3)+cos(0.52*pi*n3);y3=[x(1:1:10) zeros(1,90)];
xa (t ) [1 cos(2 100t )]cos(2 600t ) cos(2 600t ) 0.5[cos(2 700t ) cos(2 500t )]
92
• • • • • •
%cosfft0 .m N=512; n=[0:1:N-1]; fs=3000 %x=cos(2*pi/fs*100*n); x=cos(2*pi/fs*600*n)+0.5*cos(2*pi/fs*700 *n)+0.5*cos(2*pi/fs*500*n); • subplot(3,1,1);stem(n,x);ylabel('signal x(n), 0 <= n <= 511');xlabel('n') • axis([0,N-1,-2.5,2.5])
79
conv of x and h 8
7
6
5
4
y
3
2
1
0
4
6
8
10
12 14 Time index n
16
18
20
22
80
4.4 LTI系统
filter：
因果LTI系统的的零状态响应，在matlab中可以用函数
y=filter(num,den,x) 实现。 其中： num------系统传递函数分子系数组成的行向量 den-------- 系统传递函数分母系数组成的行向量 x-------输入的离散序列 y-------输出的离散序列，y的长度与x的长度一样
81
impz：
在 A(z)、B（z）已知情况下,用以下函数求系统的单 位抽样响应 h(n)， t是所记录的0到n-1的取样点数矢量： [h,t]=impz(num,den)，取样点数自动选取； [h,t]=impz(num,den, n) ，取样点数选取为0到n-1或由 矢量n 指定，取样间隔为1 ； [h,t]=impz(num,den, n，Fs) ，取样间隔为1/Fs； impz(num,den)，只在当前窗口用stem（t，h）绘出单 位抽样响应 h(n)。
•
周期信号：正弦信号，周期方波 弦信号
【例】产生一个幅度为2，频率为4Hz,相位为
A=2; f=4; phi=pi/6; w0=2*pi*f; t=0:0.01:1; x=A*sin(w0*t+phi); plot(t,x);
/6
的正
66
• 非周期信号：指数信号，阶跃信号，取样函数Sa(x)等
b=ones(1,8); n3=2;n4=9; %完整b(n)序列的非平凡值和区间端点 c=conv(a,b); nc1=n1+n3;nc2=n2+n4;%计算卷积和确定卷积非平凡区间端点 kc=nc1:nc2; %构成非平凡区间的序号自变量 stem(kc,c,'b') title('conv of x and h'); ylabel('y');xlabel('Time index n');
•%图中横坐标值即为信号的实际频率（赫兹）。 •k=0:1:N-1;w=2*pi/N*k; •subplot(3,1,3);plot(w,magX);ylabel('DFT Magnitude');xlabel('frequency in rad') •axis([0,pi,0,250]) •%图中横坐标值为信号的数字频率（弧度），乘以fs除 94 以2*pi即为信号的实际频率（赫兹）。
93
•X=fft(x,N);magX=abs(X(1:1:N)) •k=0:1:N-1;f=fs/N*k; •axis([0,1500,0,250])
•subplot(3,1,2);plot(f,magX);ylabel('DFT Magnitude');xlabel('frequency in Hz')
86
figure(1) subplot(2,1,1); stem(n,h1,'fill'); grid on; title('单位冲激响应1'); ylabel('h1(n)'); xlabel('n')
subplot(2,1,2); [h2,t]=impz(num,den, n,Fs); %计算单位冲激响应2
83
zplane：
离散系统的极－零图，在matlab中可以用以下函数作出： zplane(z,p)，已知极点、零点； zplane(num,den)； [hz，hp，ht]= zplane(z,p)，返回极点、零点、坐标轴 等的句柄。
84
【例】 已知一个系统的传递函数为
H (e
jw
0.008 0.033 e jw 0.05e j 2 w 0.033 e j 3 w 0.008e j 4 w ) 1 2.37e jw 2.7e j 2 w 1.6e j 3 w 0.41e j 4 w
96
% spectrum_com1.m &part 1 % Spectrum based on the first 10 samples of x(n)
figure(1)
n1=[0:1:9];x=cos(0.48*pi*n1)+cos(0.52*pi*n1);y1=x(1:1:10); subplot(2,1,1);stem(n1,y1);title('signal x(n), 0 <= n <= 9');xlabel('n') axis([0,10,-2.5,2.5]) Y1=fft(y1);magY1=abs(Y1(1:1:6)); k1=0:1:5;w1=2*pi/10*k1; subplot(2,1,2);stem(w1/pi,magY1);title('Samples of DTFT Magnitude'); xlabel('frequency in pi units') axis([0,1,0,10])
【例】产生一个高度为1，宽度为3，延时为2秒的矩形
脉冲 信号
t=-2:0.02:6; plot(t,rectpuls(t-2,3)); axis([-2,6,0,1.5]);
67
【例】 取样函数：用sinc(x)命令
N=1000; t=-10:20/N:10; x=sinc(t/pi); plot(t,x);grid on
signal x(n), 0 <= n <= 511
2 0 -2 0 50 100 150 200 250 n 300 350 400 450 500
DFT Magnitude
200 100 0
0
500 frequency in Hz
1000
1500
DFT Magnitude
200 100 0
0
0.5
绘出极－零图，求系统单位冲激响应h[n]，以及h[n] 的幅频相频响应图。
85
% sys_filter.m N=64; n=1:N; Fs=1024; x=[1 zeros(1,N-1)]; %产生单位冲激序列 num=[0.008 -0.033 0.05 -0.033 0.008]; den=[1 2.37 2.7 1.6 0.41]; h1=filter(num,den,x); %计算单位冲激响应1