数字信号处理B实验指导书汇总Word格式.docx

数字信号处理B实验指导书汇总Word格式.docx

《数字信号处理B实验指导书汇总Word格式.docx》由会员分享，可在线阅读，更多相关《数字信号处理B实验指导书汇总Word格式.docx（13页珍藏版）》请在冰豆网上搜索。

生成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，你看到了什么，做怎样的解释？

注意安排你的时域、频域的显示。

MATLAB程序如下

fs=[32000,24000,19000,18000,17000,16000];

m=length（fs）;

a=2*pi*6500;

b=2*pi*7000;

c=2*pi*9000;

x=1:

16;

w=（x-1）*2*pi/length（x）;

w2=0:

2*pi/256:

255*2*pi/256;

fori=1:

m

t=x/fs（i）;

y=cos（a*t）+cos（b*t）+cos（c*t）;

subplot（6,2,i*2-1）,plot（w,abs（fft（y）））;

y=[yzeros（1,256-16）];

subplot（6,2,i*2）,plot（w2,abs（fft（y）））;

共画了6*2个图,每行表示不同的采样率,每行左边表N=16的FFT,右边表示M=256的FFT

结果图如下

对于结果的讨论:

对于序列补0后的FFT与补0前相比,感觉上频谱更"

细腻"

了点,但实际上,补0后的FFT只是补0前FFT在频域上的插值而已,不会有任何更多的细节出来

另外,当fs=17000,16000时,频谱已不能分辨cos（ct）因为此时的采样频率fs小于2倍cos（c*t）的频率,这与Nyquist采样定理相符合

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等参阅联机帮助。

此题要自己编写DFT的程序

DFT的函数如下

functionout=dft（x）

m=length（x）;

out=zeros（1,m）;

t=0:

m-1;

t=t*2*pi/m;

forn=1:

out（n）=sum（x.*exp（（n-1）*i*t））;

而主程序段如下

Nmax=256;

ctime=zeros（2,Nmax-63）;

n=1;

forn=64:

tic;

ctime（1,n-63）=toc;

dft（x）;

ctime（2,n-63）=toc;

n=[64:

plot（n,ctime（1,:

）,'

b.'

holdon,plot（n,ctime（2,:

）,'

r-'

legend（'

fft'

'

dft'

FFTVSDFT'

结果

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