数字信号处理西电.docx
《数字信号处理西电.docx》由会员分享，可在线阅读，更多相关《数字信号处理西电.docx（48页珍藏版）》请在冰豆网上搜索。
数字信号处理西电
数字信号处理西电
数字信号处理上机第一次实验
实验一：
设给定模拟信号
，
的单位是ms。
（1）利用MATLAB绘制出其时域波形和频谱图（傅里叶变换），估计其等效带宽（忽略谱分量降低到峰值的3%以下的频谱）。
（2）用两个不同的采样频率对给定的
进行采样。
。
。
比较两种采样率下的信号频谱，并解释。
实验一MATLAB程序：
（1）
clc;
fs=5000;
ts=1/fs;
N=1000;
t=（-N:
N）*ts;
s=exp（-abs（t））;
plot（t,s,'linewidth',1.5）
xlabel（'时间'）
ylabel（'幅度'）
set（gca,'fontweight','b','fontsize',12）
SPL=N*100;
figure
sp=fftshift（fft（s,SPL））;
sp=sp/max（sp）*100;
freqb=-fs/2:
fs/SPL:
fs/2-fs/SPL;
plot（freqb,abs（sp））
xlabel（'频率'）
ylabel（'频谱幅度'）
set（gca,'fontweight','b','fontsize',12）
yy=abs（abs（sp）-3）;
[aa,freqind]=min（yy）;
（freqind-SPL/2）*fs/SPL
clc;
fs=1000;
ts=1/fs;
N=1000;
t=（-N:
N）*ts;
s=exp（-abs（t））;
plot（t,s,'linewidth',1.5）
xlabel（'时间'）
ylabel（'幅度'）
set（gca,'fontweight','b','fontsize',12）
SPL=N*100;
figure
sp=fftshift（fft（s,SPL））;
sp=sp/max（sp）*100;
freqb=-fs/2:
fs/SPL:
fs/2-fs/SPL;
plot（freqb,abs（sp））
xlabel（'频率'）
ylabel（'频谱幅度'）
set（gca,'fontweight','b','fontsize',12）
yy=abs（abs（sp）-3）;
[aa,freqind]=min（yy）;
（freqind-SPL/2）*fs/SPL
实验三：
设
，
，编写MATLAB程序，计算：
（1）５点圆周卷积
；
（2）６点圆周卷积
；
（3）线性卷积
；
（4）画出的
，
和
时间轴对齐。
a=[1,2,2];
b=[1,2,3,4];
y1=cconv（a,b,5）
y2=cconv（a,b,6）
y3=conv（a,b）
figure
（1）;
subplot（311）
stem（y1）;
gridon
title（'五点圆周卷积y1（n）'）;
xlabel（'n'）,ylabel（'y1（n）'）;axis（[06015]）
subplot（312）
stem（y2）;
gridon
title（'六点圆周卷积y2（n）'）;
xlabel（'n'）,ylabel（'y2（n）'）;axis（[06015]）
subplot（313）
stem（y3）;
gridon
title（'线性卷积y3（n）'）;
xlabel（'n'）,ylabel（'y3（n）'）;axis（[06015]）
实验四：
给定因果系统：
（1）求系统函数
并画出零极点示意图。
（2）画出系统的幅频特性
和相频特性
。
（3）求脉冲响应
并画序列图。
提示：
在ＭＡＴＬＡＢ中，zplane（b,a）函数可画零极点图；Freqz（b,a,N）可给出
范围内均匀间隔的
点频率响应的复振幅；Impz（b,a,N）可求
的逆变换（即脉冲响应）。
clc
a=[1,0]
b=[1,-0.9]
figure
（1）
zplane（b,a）;
title（'零极点分布图'）
w=[-3*pi:
0.01:
3*pi];
[h,phi]=freqz（b,a,w）;
figure
（2）;
subplot（2,1,1）;
plot（w,abs（h））;
gridon;
title（'幅频特性'）;
xlabel（'f/Hz'）,ylabel（'H（f）'）;
subplot（2,1,2）;
plot（w,phi）;
gridon;
title（'相频特性'）;
xlabel（'f/Hz'）,ylabel（'W（f）'）;
数字信号处理第二次实验
1.给定模拟信号
，对其进行采样，用DFT（FFT）进行信号频谱分析。
（1）确定最小采样频率和最小采样点数。
（2）若以
秒进行采样，至少需要取多少采样点？
（3）用DFT的点数
画出信号的
点DFT的幅度谱，讨论幅度谱结果。
（4）
分别为
和
，能否分辨出信号的所有频率分量。
（5）在（３）和（４）的条件下做补0FFT，分析结果。
（6）在不满足最小采样点数的情况下做补0DFT，观察是否可以分辨出两个频率分量。
程序如下：
clear
closeall
clc
%
（1）确定最小采样频率和最小采样点数
w1=4*pi;
w2=8*pi;
f1=w1/（2*pi）;
f2=w2/（2*pi）;
disp（'最小采样频率:
'）
fs1=2*max（f1,f2）;
disp（fs1）;
f=f2-f1;
disp（'最小采样点数:
'）
N=ceil（fs1/f）;
disp（N）;
%
（2）t=0.01ns采样
T=0.01;
fs2=1/T;
disp（'以t=0.01ns采样,最少采样点数为：
'）
N0=fs2/f;
disp（N0）;
%（3）（4）N=50,100,64,60时的幅度谱
w1=4*pi;
w2=8*pi;
f1=w1/（2*pi）;
f2=w2/（2*pi）;
N1=50;
N2=100;
N3=64;
N4=60;
n1=0:
N1-1;
n2=0:
N2-1;
n3=0:
N3-1;
n4=0:
N4-1;
x1=2*cos（w1*n1*T）+5*cos（w2*n1*T）;
x2=2*cos（w1*n2*T）+5*cos（w2*n2*T）;
x3=2*cos（w1*n3*T）+5*cos（w2*n3*T）;
x4=2*cos（w1*n4*T）+5*cos（w2*n4*T）;
X1=abs（fft（x1,N1））;
X2=abs（fft（x2,N2））;
X3=abs（fft（x3,N3））;
X4=abs（fft（x4,N4））;
figure
（1）
subplot（2,2,1）;
stem（n1,X1,'.'）;
title（'N=50幅度谱'）
xlabel（'n'）
ylabel（'X1'）
subplot（2,2,2）;
stem（n2,X2,'.'）;
title（'N=100幅度谱'）
xlabel（'n'）
ylabel（'X2'）
subplot（2,2,3）;
stem（n3,X3,'.'）;
title（'N=64幅度谱'）
xlabel（'n'）
ylabel（'X3'）
subplot（2,2,4）;
stem（n4,X4,'.'）;
title（'N=60幅度谱'）
xlabel（'n'）
ylabel（'X4'）
%（5）补0DFT
N5=200;
n5=0:
N5-1;
X5=abs（fft（x1,N5））;
X6=abs（fft（x2,N5））;
X7=abs（fft（x3,N5））;
X8=abs（fft（x4,N5））;
figure
（2）
subplot（2,2,1）;
stem（n5,X5,'.'）;
title（'补0后N=50幅度谱'）
xlabel（'n'）
ylabel（'X5'）
subplot（2,2,2）;
stem（n5,X6,'.'）;
title（'补0后N=100幅度谱'）
xlabel（'n'）
ylabel（'X6'）
subplot（2,2,3）;
stem（n5,X7,'.'）;
title（'补0后N=64幅度谱'）
xlabel（'n'）
ylabel（'X7'）
subplot（2,2,4）;
stem（n5,X8,'.'）;
title（'补0后N=60幅度谱'）
xlabel（'n'）
ylabel（'X8'）
%（6）N=2时不满足最小采样点数
N6=2;
n9=0:
N6-1;
x9=2*cos（w1*n9*T）+5*cos（w2*n9*T）;
X9=abs（fft（x9,N5））;
figure（3）
stem（n5,X9,'.'）;
xlabel（'n'）
ylabel（'X9'）
title（'N=2时补0后的幅度谱'）
运行结果：
2.设雷达发射线性调频信号
，
，采样率
，采样点数
。
回波信号
，
，
。
（1）画出
的频谱。
（2）利用DFT的时延性质产生
，比较直接在时域产生和在频域产生（再变换到时域）的结果是否相同。
（3）匹配滤波的结果是
，（“
”表示线性卷积）。
分别用直接线性卷积和DFT的卷积定理求解
。
比较二者结果，并记录两种方法的运行时间（用tic，toc指令）。
（4）画出
的频谱。
程序如下：
figure;
plot（[-0.5:
1/（fft_num）:
0.5-1/（fft_num）],...
fftshift（20*log10（abs（fft（ht,fft_num）））））%将线性调频信号转换到频域并将零频搬至频谱中央
axis（[-0.50.51050]）
xlabel（'归一化频率'）
ylabel（'幅度/dB'）
title（'h（t）频谱'）
第二问
时域构造回波信号，时延通过补零实现
s_shiyu=[zeros（1,shiyan1*fs）,ht,zeros（1,N-shiyan1*fs）]+[zeros（1,shiyan2*fs）,ht,zeros（1,N-shiyan2*fs）];
figure;
plot（[0:
2*N-1],abs（s_shiyu））
axis（[02*N-102.5]）
xlabel（'距离单元'）
ylabel（'幅度'）
title（'时域法s（t）'）
figure;
plot（[-0.5:
1/（fft_num）:
0.5-1/（fft_num）],...
fftshift（20*log10（abs（fft（s_shiyu,fft_num）））））%将回波信号转换到频域并将零频搬至频谱中央
axis（[-0.50.51060]）
xlabel（'归一化频率'）
ylabel（'幅度/dB'）
title（'时域法s（t）频域'）
频域构造回波信号，时延通过DFT时延性质产生
L=2*N;
P=fft（ht,L）;
P_shiyan1=P.*exp（-j*2*pi*fs*[0:
L-1]/L*shiyan1）;%目标1频谱
P_shiyan2=P.*exp（-j*2*pi*fs*[0:
L-1]/L*shiyan2）;%目标2频谱
s_pinyu=ifft（P_shiyan1）+ifft（P_shiyan2）;
figure;
plot（[0:
2*N-1],abs（s_pinyu））
axis（[02*N-102.5]）
xlabel（'距离单元'）
ylabel（'幅度'）
title（'频域法s（t）'）
figure;
plot（[-0.5:
1/（fft_num）:
0.5-1/（fft_num）],...
fftshift（20*log10（abs（fft（s_pinyu,fft_num）））））
axis（[-0.50.51050]）
axis（[-0.50.51060]）
xlabel（'归一化频率'）
ylabel（'幅度/dB'）
title（'频域法s（t）频域'）
figure;
plot（[0:
2*N-1],abs（s_shiyu-s_pinyu））
axis（[02*N-101]）
xlabel（'距离单元'）
ylabel（'幅度'）
title（'时域法与频域法回波之差'）
第三问
tic
y_shiyu=conv（s_pinyu,conj（fliplr（ht）））;%%%时域匹配滤波
toc
y_shiyu_quchu=y_shiyu（1,N:
end）;%%%%去暂态点
figure;
plot（[0:
2*N-1],abs（y_shiyu_quchu））
xlabel（'距离单元'）
ylabel（'幅度'）
title（'时域匹配滤波'）
tic
y_pinyu=ifft（fft（s_pinyu）.*conj（P））;频域匹配滤波
toc
figure;
plot（[0:
2*N-1],abs（y_pinyu））
xlabel（'距离单元'）
ylabel（'幅度'）
title（'频域匹配滤波'）
第四问
figure
plot（[-0.5:
1/L:
0.5-1/L],fftshift（20*log10（abs（fft（y_pinyu）））））
xlabel（'归一化频率'）
ylabel（'幅度/dB'）
title（'y（t）频域'）运行结果：
Elapsedtimeis0.235671seconds.
Elapsedtimeis0.001764seconds.
>>
数字信号处理上机第三次实验
1.IR滤波器设计
（1）用matlab确定一个数字IIR低通滤波器所有四种类型的最低阶数。
指标如下：
40kHz的采样率，4kHz的通带边界频率，8kHz的阻带边界频率，0.5dB的通带波纹，40dB的最小阻带衰减。
并在同一张图中画出每种w。
（2）用matlab确定一个数字IIR高通滤波器所有四种类型的最低阶数。
指标如下：
3500Hz的采样率，1050Hz的通带边界频率，600Hz的阻带边界频率，1dB的通带波纹，50dB的最小阻带衰减。
并在同一张图中画出每种滤波器的频率响应。
（3）用matlab确定一个数字IIR带通滤波器所有四种类型的最低阶数。
指标如下：
7kHz的采样率，1.4kHz和2.1kHz的通带边界频率，1.05kHz和2.45kHz的阻带边界频率，0.4dB的通带波纹，50dB的最小阻带衰减。
并在同一张图中画出每种滤波器的频率响应。
（4）用matlab确定一个数字IIR带阻滤波器所有四种类型的最低阶数。
指标如下：
12kHz的采样率，2.1kHz和4.5kHz的通带边界频率，2.7kHz和3.9kHz的阻带边界频率，0.6dB的通带波纹，45dB的最小阻带衰减。
并在同一张图中画出每种滤波器的频率响应。
用到的函数：
butter，buttord，cheb2ord，chebl1，cheby2，
ellip，ellipord.
程序如下：
（1）clc
clearall
closeall
fc=40;fp=4;fs=8;rp=0.5;rs=40;
wp=2*pi*fp/fc;
ws=2*pi*fs/fc;
disp（'Inbuttord'）
[n,wc]=buttord（wp,ws,rp,rs,'s'）
[b,a]=butter（n,wc,'low','s'）;
w=0:
0.001:
6;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
plot（w,h,'r-'）
disp（'Incheb1ord'）
[n,wpo]=cheb1ord（wp,ws,rp,rs,'s'）
[b,a]=cheby1（n,rp,wpo,'low','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'b-'）
disp（'Incheb2ord'）
[n,wso]=cheb2ord（wp,ws,rp,rs,'s'）
[b,a]=cheby2（n,rs,wso,'low','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'k-'）
disp（'Inellipord'）
[n,wc]=ellipord（wp,ws,rp,rs,'s'）
[b,a]=ellip（n,rp,rs,wc,'low','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'m-'）
title（'滤波器的频率响应'）
legend（'巴特沃斯','切比雪夫Ⅰ型','切比雪夫Ⅱ型','椭圆'）
gridon
xlabel（'w'）
ylabel（'h'）
运行结果：
巴特沃斯
n=
9
wc=
0.7533
切比雪夫Ⅰ型
n=
5
wpo=
0.6283
切比雪夫Ⅱ型
n=
5
wso=
1.2069
椭圆
n=
4
wc=
0.6283
（2）
clc
clearall
closeall
fc=3500;fp=1050;fs=600;rp=1;rs=50;
wp=2*pi*fp/fc;
ws=2*pi*fs/fc;
disp（'Inbuttord'）
[n,wc]=buttord（wp,ws,rp,rs,'s'）
[b,a]=butter（n,wc,'high','s'）;
w=0:
0.001:
6;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
plot（w,h,'r-'）
disp（'Incheb1ord'）
[n,wpo]=cheb1ord（wp,ws,rp,rs,'s'）
[b,a]=cheby1（n,rp,wpo,'high','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'b-'）
disp（'Incheb2ord'）
[n,wso]=cheb2ord（wp,ws,rp,rs,'s'）
[b,a]=cheby2（n,rs,wso,'high','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'k-'）
disp（'Inellipord'）
[n,wc]=ellipord（wp,ws,rp,rs,'s'）
[b,a]=ellip（n,rp,rs,wc,'high','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'m-'）
title（'滤波器的频率响应'）
legend（'巴特沃斯','切比雪夫Ⅰ型','切比雪夫Ⅱ型','椭圆'）
axis（[0,6,-100,0]）
gridon
xlabel（'w'）
ylabel（'h'）
运行结果：
巴特沃斯
n=
12
wc=
1.7402
切比雪夫Ⅰ型
n=
7
wpo=
1.8850
切比雪夫Ⅱ型
n=
7
wso=
1.2049
椭圆
n=
5
wc=
1.8850
（3）
clc
clearall
closeall
fc=7;fp=[1.4,2.1];fs=[1.05,2.45];rp=0.4;rs=50;
wp=2*pi*fp/fc;
ws=2*pi*fs/fc;
disp（'Inbuttord'）
[n,wc]=buttord（wp,ws,rp,rs,'s'）
[b,a]=butter（n,wc,'s'）;
w=0:
0.001:
6;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
plot（w,h,'r-'）
disp（'Incheb1ord'）
[n,wpo]=cheb1ord（wp,ws,rp,rs,'s'）
[b,a]=cheby1（n,rp,wpo,'s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'b-'）
disp（'Incheb2ord'）
[n,wso]=cheb2ord（wp,ws,rp,rs,'s'）
[b,a]=cheby2（n,rs,wso,'s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'k-'）
disp（'Inellipord'）
[n,wc]=ellipord（wp,ws,rp,rs,'s'）
[b,a]=ellip（n,rp,rs,wc,'s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
holdon
plot（w,h,'m-'）
title（'滤波器的频率响应'）
legend（'巴特沃斯','切比雪夫Ⅰ型','切比雪夫Ⅱ型','椭圆'）
axis（[0,6,-100,20]）
gridon
xlabel（'w'）
ylabel（'h'）
运行结果：
巴特沃斯
n=
12
wc=
1.23051.9250
切比雪夫Ⅰ型
n=
7
wpo=
1.25661.8850
切比雪夫Ⅱ型
n=
7
wso=
1.10502.1437
椭圆
n=
5
wc=
1.25661.8850
（4）
clc
clearall
closeall
fc=12;fp=[2.1,4.5];fs=[2.7,3.9];rp=0.5;rs=40;
wp=2*pi*fp/fc;
ws=2*pi*fs/fc;
disp（'Inbuttord'）
[n,wc]=buttord（wp,ws,rp,rs,'s'）
[b,a]=butter（n,wc,'stop','s'）;
w=0:
0.001:
6;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
plot（w,h,'r-'）
disp（'Incheb1ord'）
[n,wpo]=cheb1ord（wp,ws,rp,rs,'s'）
[b,a]=cheby1（n,rp,wpo,'stop','s'）;
[h,w]=freqs（b,a,w）;
h=20*log10（abs（h））;
升级会员 