数字信号处理的MATLAB实现

昆明理工大学信息工程与自动化学院学生实验报告

（2011—2012 学年第二学期）

课程名称：数字信号处理开课实验室：信自楼111 2012 年 5 月 31 日年级、专业、班生医学号姓

名

成绩

实验项目名称数字信号处理的matlab 实现指导教师

教

师

评语教师签名：

年月日

一．实验目的

熟练掌握matlab的基本操作。

了解数字信号处理的MATLAB实现。

二．实验设备

安装有matlab的PC机一台。

三．实验内容

.1.求信号x(n)=cos(6.3Пn/3)+cos(9.7Пn/30)+cos(15.3Пn/30),0≤n≤29的幅度频谱.

2. 用冲击响应不变法设计一个Butterworth低通数字滤波器,要求参数为:

Wp=0.2Пαp=1dB

Ws=0.3Пαs=15dB

3.用双线性变换法设计一个Chebyshev高通IIR滤波器,要求参数为:

Wp=0.6Пαp=1dB

Ws=0.4586Пαs=15dB

4.用窗函数法设计一个低通FIR滤波器,要求参数为:

Wp=0.2Пαp=0.3dB

Ws=0.25Пαs=50dB

5.用频率抽样法设计一个带通FIR滤波器,要求参数为:

W1s=0.2П W1p=0.35П W2p=0.65П W2s=0.8П

αs=60dB αp=1dB

6.根据 4 点矩形序列，( n ) = [1 1 1 1] 。做 DTFT 变换，再做 4 点 DFT 变换。然后分别补零做 8 点 DFT 及 16 点 DFT。

7.调用filter解差分方程，由系统对u(n)的响应判断稳定性

8编制程序求解下列系统的单位冲激响应和阶跃响应。 y[n]+ 0.75y[n -1]+ 0.125y[n -2] = x[n]- x[n -1]

四．实验源程序

1.

n=[0:1:29];

x=cos(6.3*pi*n/30)+cos(9.7*pi*n/30)+cos(15.3*pi*n/30);

subplot(2,2,1);

stem(n,x);

title('信号x(n),0<=n<=29');

X=fft(x);

magX=abs(X(1:1:16));

k=0:1:15;

subplot(2,2,2);

stem(k,magX);

title('FFT幅度');

2.

wp=0.2*pi; %通带边界数字频率（Hz）

ws=0.3*pi; %阻带边界数字频率（Hz）

Rp=1; %通带波动（dB）

As=15; %阻带波动（dB）

T=1; %取T=1

OmegaP=wp*T; %模拟原型通带边界频率

OmegaS=ws*T; %模拟原型阻带边界频率

[cs,ds]=afd_butt(OmegaP,OmegaS,Rp,As);%调用模拟Butterworth滤波器的设计

[b,a]=imp_invr(cs,ds,T); %调用模拟到数字滤波器转换函数%[b0,B,A]=dir2cas(b,a); %直接型转换为级联型系数的计算[H,w]=freqz(b,a,501); %调用计算频率特性函数

mag=abs(H);pha=angle(H);

subplot(2,2,1);plot(w/pi,mag);title('幅频特性');

subplot(2,2,2);plot(w/pi,-pha/pi);title('相频特性');

3.

ws=0.4586*pi; %通带边界数字频率（Hz)

wp=0.6*pi; %阻带边界数字频率（Hz）

Rp=1; %通带波动（dB）

As=15; %阻带衰减（dB）

[N,wn]=cheb1ord(wp/pi,ws/pi,Rp,As);%Chebyshew滤波器参数计算

[b,a]=cheby1(N,Rp,wn,'high'); %数字Chebyshew高通滤波器设计

%[b0,B,A]=dir2cas(b,a); %直接型转换为级联型系数计算

[H,w]=freqz(b,a,501); %调用计算频率特性函数

mag=abs(H);pha=angle(H); %计算幅频和相频特性

subplot(2,2,1);plot(w/pi,mag);title('幅频特性');

subplot(2,2,2);plot(w/pi,-pha/pi);title('相频特性');

4.

p=0.2*pi;ws=0.3*pi;

tr_width=ws-wp;

M=ceil(6.6*pi/tr_width)+1;

n=[0:1:M-1];

wc=(ws+wp)/2;

hd=ideal_lp(wc,M);

w_ham=(hamming(M))';

h=hd.*w_ham;

[H,w]=freqz(h,[1],501);

mag=abs(H);pha=angle(H);

db=20*log10((mag+eps)/max(mag));

delta_w=2*pi/1000;

Rp=-(min(db(1:1:wp/delta_w+1)));

As=-round(max(db(ws/delta_w+1:1:501)));

subplot(2,2,1);stem(n,h);title('单位抽样响应');

subplot(2,2,2);plot(w/pi,db);title('幅频特性(db)');

Rp=-(min(db(1:1:wp/delta_w+1)))

As=-round(max(db(ws/delta_w+1:1:501)))

5.

M=40;alpha=(M-1)/2;l=0:M-1;wl=(2*pi/M)*l;

T1=0.109021;T2=0.59417456;

Hrs=[zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1, 4)];

Hdr=[0,0,1,1,0,0];wdl=[0,0.2,0.35,0.65,0.8,1];

k1=0:floor((M-1)/2);k2=floor((M-1)/2)+1:M-1;

angH=[-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2)];

H=Hrs.*exp(j*angH);

h=real(ifft(H,M));

[db,mag,pha,grd,w]=freqz_m(h,1);

[Hr,ww,a,L]=Hr_Type2(h);

subplot(2,2,1);stem(l,h);title('单位抽样响应');

subplot(2,2,2);plot(w/pi,db);axis([0,1,-100,10]);

title('幅频特性')

set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1]);

6.

% x(n)=[1 1 1 1],长度为4的离散序列,求DTFT[x(n)];

n=0:3;

x=[1 1 1 1];

k=0:1000;

w=(2*pi/1000)*k;

X=x*(exp(-j*2*pi/1000)).^(n'*k);

magX=abs(X);angX=angle(X);

realX=real(X);imagX=imag(X);

hold on

figure(1);

subplot(2,1,1);plot(k/500,magX);grid