dsp实验数字信号处理丁玉美第三版精.docx

dsp实验数字信号处理丁玉美第三版精.docx

《dsp实验数字信号处理丁玉美第三版精.docx》由会员分享，可在线阅读，更多相关《dsp实验数字信号处理丁玉美第三版精.docx（13页珍藏版）》请在冰豆网上搜索。

dsp实验数字信号处理丁玉美第三版精

第１章离散时间信号与系统

例1-1用MATLAB计算序列{-201–13}和序列{120-1}的离散卷积。

解MATLAB程序如下：

a=[-201-13];

b=[120-1];

c=conv（a,b）;

M=length（c）-1;

n=0:

1:

M;

stem（n,c）;

              xlabel（'n'）;ylabel（'幅度'）;

图1.1给出了卷积结果的图形，求得的结果存放在数组c中为：

{-2-413151-3}。

例1-2用MATLAB计算差分方程

当输入序列为时的输出结果。

解MATLAB程序如下：

N=41;

a=[0.8-0.440.360.22];

b=[10.7-0.45-0.6];

x=[1zeros（1,N-1）];

k=0:

1:

N-1;

y=filter（a,b,x）;

stem（k,y）

xlabel（'n'）;ylabel（'幅度'）

          图1.2给出了该差分方程的前41个样点的输出，即该系统的单位脉冲响应。

例1-3用MATLAB计算例1-2差分方程

所对应的系统函数的DTFT。

解例1-2差分方程所对应的系统函数为：

       其DTFT为

     用MATLAB计算的程序如下：

k=256;

num=[0.8-0.440.360.02];

den=[10.7-0.45-0.6];

w=0:

pi/k:

pi;

h=freqz（num,den,w）;

subplot（2,2,1）;

plot（w/pi,real（h））;grid

title（'实部'）

xlabel（'\omega/\pi'）;ylabel（'幅度'）

subplot（2,2,2）;

plot（w/pi,imag（h））;grid

title（'虚部'）

xlabel（'\omega/\pi'）;ylabel（'Amplitude'）

subplot（2,2,3）;

plot（w/pi,abs（h））;grid

title（'幅度谱'）

xlabel（'\omega/\pi'）;ylabel（'幅值'）

subplot（2,2,4）;

plot（w/pi,angle（h））;grid

title（'相位谱'）

xlabel（'\omega/\pi'）;ylabel（'弧度'）

第2章离散傅里叶变换及其快速算法 

例2-1对连续的单一频率周期信号按采样频率采样，截取长度N分别选N=20和N=16，观察其DFT结果的幅度谱。

解此时离散序列，即k=8。

用MATLAB计算并作图，函数fft用于计算离散傅里叶变换DFT，程序如下：

k=8;

n1=[0:

1:

19];

xa1=sin（2*pi*n1/k）;

subplot（2,2,1）

plot（n1,xa1）

xlabel（'t/T'）;ylabel（'x（n）'）;

xk1=fft（xa1）;xk1=abs（xk1）;

subplot（2,2,2）

stem（n1,xk1）

xlabel（'k'）;ylabel（'X（k）'）;

n2=[0:

1:

15];

xa2=sin（2*pi*n2/k）;

subplot（2,2,3）

plot（n2,xa2）

xlabel（'t/T'）;ylabel（'x（n）'）;

xk2=fft（xa2）;xk2=abs（xk2）;

                       subplot（2,2,4）

stem（n2,xk2）

xlabel（'k'）;ylabel（'X（k）'）;

        计算结果示于图2.1，（a）和（b）分别是N=20时的截取信号和DFT结果，由于截取了两个半周期，频谱出现泄漏；（c）和（d）分别是N=16时的截取信号和DFT结果，由于截取了两个整周期，得到单一谱线的频谱。

上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。

例2-2用FFT计算两个序列

的互相关函数。

解用MATLAB计算程序如下：

x=[13-112331];

y=[21-1120-13];

k=length（x）;

xk=fft（x,2*k）;

yk=fft（y,2*k）;

rm=real（ifft（conj（xk）.*yk））;

rm=[rm（k+2:

2*k）rm（1:

k）];

m=（-k+1）:

（k-1）;

stem（m,rm）

xlabel（'m'）;ylabel（'幅度'）;

                   其计算结果如图2.2所示。

例2-3计算两个序列的的互相关函数，其中

x（n）={23521–100123530–1–2012}

y（n）=x（n-4）+e（n）,e（n）为一随机噪声，在MATLAB中可以用随机函数rand产生

解用MATLAB计算程序如下：

x=[23521-100123530-1-2012];

y=[000023521-100123530-1-2012];

k=length（y）;

e=rand（1,k）-0.5;

y=y+e;

xk=fft（x,2*k）;

yk=fft（y,2*k）;

rm=real（ifft（conj（xk）.*yk））;

rm=[rm（k+2:

2*k）rm（1:

k）];

m=（-k+1）:

（k-1）;

stem（m,rm）

xlabel（'m'）;ylabel（'幅度'）;

计算结果如图2.3（a），我们看到最大值出现在m=4处，正好是y（n）对于x（n）的延迟。

2.3（b）是x（n）的自相关函数，他和y（n）的区别除时间位置外，形状也略不同，这是由于y（n）受到噪声的干扰。

第3章无限长单位脉冲响应（IIR）滤波器的设计方法

例3-1设采样周期T=250μs（采样频率fs=4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc=1kHz。

[B,A]=butter（3,2*pi*1000,'s'）;

[num1,den1]=impinvar（B,A,4000）;

[h1,w]=freqz（num1,den1）;

[B,A]=butter（3,2/0.00025,'s'）;

[num2,den2]=bilinear（B,A,4000）;

[h2,w]=freqz（num2,den2）;

f=w/pi*2000;

plot（f,abs（h1）,'-.',f,abs（h2）,'-'）;

grid;

xlabel（'频率/Hz'）

ylabel（'幅值/dB'）

  程序中第一个butter的边界频率2π×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图3.1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。

脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。

同时，也看到双线性变换法，在z=-1即ω=π或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在Ω=∞处的三阶传输零点通过映射形成的。

例3-2设计一数字高通滤波器，它的通带为400～500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1,000Hz。

wc=2*1000*tan（2*pi*400/（2*1000））;

wt=2*1000*tan（2*pi*317/（2*1000））;

[N,wn]=cheb1ord（wc,wt,0.5,19,'s'）;

[B,A]=cheby1（N,0.5,wn,'high','s'）;

[num,den]=bilinear（B,A,1000）;

[h,w]=freqz（num,den）;

f=w/pi*500;

plot（f,20*log10（abs（h）））;

axis（[0,500,-80,10]）;

grid;

xlabel（''）

ylabel（'幅度/dB'）

图3.2给出了MATLAB计算的结果，可以看到模拟滤波器在Ω=∞处的三阶零点通过高通变换后出现在ω=0（z=1）处，这正是高通滤波器所希望得到的。

例3-3设计一巴特沃兹带通滤波器，其３dB边界频率分别为f2=110kHz和f1=90kHz，在阻带f3=120kHz处的最小衰减大于１０dB，采样频率fs=400kHz。

w1=2*400*tan（2*pi*90/（2*400））;

w2=2*400*tan（2*pi*110/（2*400））;

wr=2*400*tan（2*pi*120/（2*400））;

[N,wn]=buttord（[w1w2],[0wr],3,10,'s'）;

[B,A]=butter（N,wn,'s'）;

[num,den]=bilinear（B,A,400）;

[h,w]=freqz（num,den）;

f=w/pi*200;

plot（f,20*log10（abs（h）））;

axis（[40,160,-30,10]）;

grid;

xlabel（'频率/kHz'）

ylabel（'幅度/dB'）

图3.3给出了MATLAB计算的结果，可以看出数字滤波器将无穷远点的二阶零点映射为z=±1的二阶零点，数字带通滤波器的极点数是模拟低通滤波器的极点数的两倍。

例3-4一数字滤波器采样频率fs=1kHz，要求滤除100Hz的干扰，其３dB的边界频率为95Hz和105Hz，原型归一化低通滤波器为

w1=95/500;

w2=105/500;

[B,A]=butter（1,[w1,w2],'stop'）;

[h,w]=freqz（B,A）;

f=w/pi*500;

plot（f,20*log10（abs（h）））;

axis（[50,150,-30,10]）;

grid;

xlabel（'频率/Hz'）

ylabel（'幅度/dB'）

           图3.4为MATLAB的计算结果

第4章有限长单位脉冲响应（FIR）滤波器的设计方法

例2用凯塞窗设计一FIR低通滤波器，低通边界频率，阻带边界频率，阻带衰减不小于