数字信号处理B实验指导书Word文档格式.docx
《数字信号处理B实验指导书Word文档格式.docx》由会员分享,可在线阅读,更多相关《数字信号处理B实验指导书Word文档格式.docx(11页珍藏版)》请在冰豆网上搜索。
n阶的均匀分布随机阵;
Y=rand(m,n)生成m×
n阶的随机阵;
rand返回在[0,1]区间上的一个随机数;
将上面的rand写成randn则可以生成均值为0、方差为1的正态分布的随机变量。
6.全1矩阵生成函数ones(m,n):
生成m×
n阶全1矩阵
7.全0矩阵生成函数zeros(m,n):
n阶全0矩阵
8.离散序列绘图函数stem
stem(y)以1、2、3…为横坐标,
y为纵坐标画杆形图;
stem(x,y)以x为横坐标,
y为纵坐标画杆形图(x与y数据个数必须一致);
stem(…,’fill’)选项’fill’指定杆顶为实心,若无此选项则默认空心。
stem(...,LineSpec)参数LineSpec指定杆形图的线形、数据标志符号及颜色,具体用法可查看MATLAB帮助
9.线性坐标平面绘图函数plot
用法与stem
类似,具体用法可查看MATLAB帮助
以上为MATLAB内置函数(在此仅为同学复习MATLAB提供)
四、在MATLAB程序中变量赋值注意问题
在MATLAB
中,对变量赋值时其维数可以按需要动态地改变,这样虽然方便程序设计但同时容易出错。
另外,频繁分配变量空间会大大降低程序的执行速度,因而应该尽量避免不必要的矩阵、向量维数的改变。
通常先用zeros()函数给变量分配足够大小的空间,再对变量进行赋值。
例:
依次执行下面的语句
tic%开始计时
fori=1:
10000
c(i)=i;
%每次都重新分配空间
end
toc%读取计时时间
d=zeros(1,10000);
%预先分配空间
d(i)=i;
%直接赋值,不必重新分配空间
运行结果如下:
elapsed_time=
1.1560
0.0470
从结果可以看出,第2种赋值方法所用的时间比第1种方法所用时间少得多(以上是在主频为2.66GHZ的机器上运行的结果)。
实验2离散信号、系统的频域表示
1.考察抽样间隔对信号频谱的影响;
2.掌握用FFT做谱分析的方法。
1、用DFT/FFT对模拟信号做傅里叶分析
以频率fs对以下信号抽样N点
xa(t)=cos(at)+cos(bt)+cos(ct)
相应的参数是
a=2*pi*6500,b=2*pi*7000,c=2*pi*9000
fs=32000,N=16
对这N点序列作N点DFT,观察其幅频特性,如果
X=fft(x)
w是频率坐标向量,你可以考虑用stem(w,abs(X)),plot(w,abs(X)),plot(w,abs(X),'
*'
)来显示,然后确定用哪种显示方式。
接下来,对x作M=256点FFT,这意味着在x后补了M-N个0,再观察幅频特性;
现在令抽样点数N=256,再对这个抽样序列作N点FFT,观察幅频特性。
通过这些观察,你能作出什么判断?
试着改变fs,令其分别为24000,19000,18000,17000,16000,你看到了什么,做怎样的解释?
注意安排你的时域、频域的显示。
2、参考以下程序段,对DFT和FFT算法做计算效率比较
Nmax=2048;
fft_time=zeros(1,Nmax);
forn=1:
Nmax
x=rand(1,n);
t=clock;
fft(x);
fft_time(n)=etime(clock,t);
n=[1:
Nmax];
plot(n,fft_time,'
.'
xlabel('
N'
);
ylabel('
TimeinSec.'
title('
FFTexecutiontimes'
比如,从N1点到N2点内的效率比较,设N1=8,N2可以考虑取256,或更大些。
大致需要这样一个循环
forn=N1:
N2
产生数据(可以用固定数据,或是随机数据)
计时开始
DFT
计时结束,统计
FFT
end
显示统计数据
计时函数clock,tic等参阅联机帮助。
3、对抽取/内插前后的信号做傅里叶分析
本次实验的信号均假定起始时间下标为0,也就是对信号x(n)作N:
1抽取,只要
y=x(1:
N:
end)
即可,而1:
N内插则为
y=zeros(1,N*length(x));
y([1:
length(y)])=x;
我们要着重观察的是抽取、内插后频谱的改变。
本实验的数据放在updn.mat文件内,执行语句
loadupdn
要用的数据就会载入数组siga和sigb,如何获取数组大小的信息?
对siga做抽取,sigb做内插实验,试用
N=2,3,4
做抽取,内插,观察它们频谱的变化。
提示:
做谱分析时补一定的0在序列尾部。
实验3FIR数字滤波器的设计
掌握用MATLAB设计FIR数字滤波器的方法,重点要掌握窗函数法。
二、实验原理
用MATLAB提供的函数,设计FIR数字滤波器
三、实验内容与要求
1.参考示范程序
窗法:
Examples7.8-11(在E:
\DSP-A\RefProgram\CHAP_07目录下)
最优化设计(Parks-McClellan-Remez):
Examples7.23-25(在E:
在阅读这些例题时,应该注意:
A)如何从滤波器指标要求导出其它设计参数,如确定窗口类型、滤波器的阶数等;
B)例题中验证设计得到的滤波器是否满足设计指标的语句;
C)特别是在最优化设计的例中的迭代设计过程;
D)把滤波器的设计和其特性的显示分开,你自己的显示不一定要照搬书上的。
特别最优化设计的那几个例子程序,显示部分需要函数ampl_res(),需要自己写。
在这基础上,试着直接用我们讲的MATLAB函数fir1()、fir2()进行设计
2.用Kaiser窗设计一个FIR数字带阻滤波器,对模拟信号
a=2*pi*6500,b=2*pi*7000,c=2*pi*9000
滤波,要求滤去7000Hz的频率成分。
系统采样率为
fs=32000Hz
这同我们第二次实验。
但采样点数应该比较大,可以用N=4096。
滤波器的Rp=0.25dB,As=50dB,过渡带宽可以用模拟频率(例如200Hz)也可以用数字频率指定。
还可以改变As(比如30dB)观察滤波效果。
这个实验一般不须在时域观察,主要应在频域作观察。
用
freqz(x,1)
就能观察有限长序列的DTFT,对对数显示不习惯的,可以用
[H,w]=freqz(x,1);
得到数据,再用
freqzplot(H,w,'
linear'
等方式显示。
具体用法参见参考程序。
3.最优化设计方法为选做项目。
四、参考程序:
1、FIR1Demo1.m用fir1函数设计高通滤波器例,各种窗口都试了一下;
2、FIR1Demo2.m用fir1函数设计多通带滤波器例;
3、FIR2Demo3.m用fir2函数设计多通带滤波器例;
4、FIRLsRemesDemo.m用firls和remez函数设计例,fir2也在这里做了对照。
1和4项所列的文件注释详细些。
FIR1Demo1.m
%用窗函数设计线性相位FIR滤波器
%教科书讲到的6种窗都用一下
%本例未涉及如何确定滤波器的设计参数
figNo=1;
%显示用窗口号,接连用了3个
n=48%滤波器阶数
Wn=0.35;
%截止频率
beta=0.1102*(60-8.7);
%Kaiser窗参数,假设阻带要求60分贝衰减(p.253)
%生成窗口矩阵,各窗函数看教科书pp.243-53
%先创建一个空矩阵,然后再把各窗函数列向量加进去
Windows=[];
Windows=[Windows,rectwin(n+1)];
Windows=[Windows,bartlett(n+1)];
Windows=[Windows,blackman(n+1)];
Windows=[Windows,hann(n+1)];
Windows=[Windows,hamming(n+1)];
Windows=[Windows,kaiser(n+1,beta)];
%接下来用一个循环完成用各窗口函数设计的滤波器
%并得到相应的传输函数向量,按列放在矩阵Hs中
%在这个循环中,win依次取矩阵Windows的各列
Hs=[];
forwin=Windows
b=fir1(n,Wn,'
high'
win);
%把high参数去掉就是低通滤波器
[h,w]=freqz(b,1);
Hs=[Hs,h];
%现在用3种尺度来显示,请观察对比各窗口
figure(figNo)
freqzplot(Hs,w,'
figure(figNo+1)
squared'
figure(figNo+2)
freqzplot(Hs,w)%用分贝数显示
figure(figNo)%把第一个窗口推到前头
FIR1Demo2.m
%用窗函数设计线性相位FIR带通/带阻滤波器
n=48;
Wn=[0.350.65];
%通带或阻带
%Hamingwindowbydefault
b=fir1(n,Wn,'
stop'
%再省去第三个参数'
就是带通
win=rectwin(n+1);
%矩形窗,宽度是滤波器阶数+1
bb=fir1(n,Wn,'
[H,w]=freqz(b,1,512);
[HH]=freqz(bb,1,512);
TH=[H,HH];
freqzplot(TH,w,'
FIR2Demo3.m
%设计2个带通/带阻滤波器演示
n=48
f=[0.20.40.60.8];
%f的元素一定不能包含0和1!
bhm=fir1(n,f,'
DC-1'
%'
DC-0'
表示频率[0,0.2]是阻带,'
是通带
win=rectwin(n+1);
%上面用缺省海明窗,下面用的矩形窗
brec=fir1(n,f,'
%显示
[Hm,w]=freqz(bhm,1);
[Hrec]=freqz(brec,1);
H2=[Hm,Hrec];
freqzplot(H2,w,'
FIRLsRemesDemo.m
%最优化FIR线性相位滤波器设计例
%均方误差最小准则
%b=firls(n,f,m)
%最大误差最小化准则,雷米兹交换算法
%b=remez(n,f,m)
%参数:
%b:
返回的滤波器传输函数分子多项式的系数,也就是单位样本响应
%n:
滤波器的阶数
%f:
行向量,以0开始递增,各关键频率点,以1结束
%m:
行向量,维数与f相同,各关键频率点对应的响应幅度值
%n=20;
%Filterorder
%f=[00.40.51];
%Frequencybandedges
%m=[1100];
%Desiredamplitudes
%下面的例子是用这两种方法设计的两个通带的滤波器
%再加上FIR2设计的作为对照
%useplot是控制显示方式的变量,参看下面显示部分
useplot=0;
%试着改动各参数观察一下,比如改变过渡带的宽度、阶数等
n=40;
m=[11001100];
f=[00.30.40.50.60.70.81];
b1=firls(n,f,m);
b2=remez(n,f,m);
%用FIR2设计
b3=fir2(n,f,m);
%下面的代码显示比较它们的频率响应
[Hb1,w]=freqz(b1,1);
%分母多项式只有常数项1,没指定返回点数,缺省512点
[Hb2]=freqz(b2);
%这里把分母项省略了,返回的点数和上行一样,所以不再需要w
Hb3=freqz(b3);
ifuseplot
%用一句plot画3条曲线,第2条用绿色,虚线,第3条红色;
grid画出网格线
subplot(2,1,1)
plot(w/pi,abs(Hb1),w/pi,abs(Hb2),'
g--'
w/pi,abs(Hb3),'
r'
),grid
subplot(2,1,2)
%plot(w/pi,angle(Hb),w/pi,angle(Hbb),'
r--'
plot(w/pi,unwrap(angle(Hb1)),w/pi,unwrap(angle(Hb2)),'
w/pi,unwrap(angle(Hb3)),'
grid
else%也可以用下面的来显示
%Hb,Hbb是和w同维数的列向量(从Wokspace里可以看到)
%H是2列合成的矩阵
H=[Hb1,Hb2,Hb3];
%用不同的方式显示,如果三种都要同时看,你还得加几句
%例如用pause在时间上隔开,或是用3个窗口
figure
(1)
freqzplot(H,w)
figure
(2)
freqzplot(H,w,'
figure(3)
五、实验所用数值处理函数表
MATLAB中提供了一些数值处理函数,在编程时常会用到,它们也是按照数组运算规则进行运算的。
函数名
功能
fix
向零取整
ceil
向正无穷方向取整
mod
除法求余
(与除数同号)
round
四舍五入
floor
向负无穷方向取整
rem
(与被除数同号)
我们的网站:
出品:
天涯工作室
联系我们:
2447417264