数字信号处理实验报告(自己的实验报告)

数字信号处理

实验报告

西南交通大学信息科学与技术学院

姓名：伍先春

学号：20092487

班级：自动化1班

指导老师：张翠芳

实验一序列的傅立叶变换

实验目的

进一步加深理解DFS,DFT 算法的原理;研究补零问题;快速傅立叶

变换（FFT ）的应用。

实验步骤

1. 复习DFS 和DFT 的定义，性质和应用；

2. 熟悉MATLAB 语言的命令窗口、编程窗口和图形窗口的使用；利用提

供的程序例子编写实验用程序；按实验内容上机实验，并进行实验结果分析；写出完整的实验报告，并将程序附在后面。

实验内容

1. 周期方波序列的频谱试画出下面四种情况下的的幅度频谱,并分析补零后，对信号频谱的影响。

2. 有限长序列x(n)的DFT

（1） 取x(n)(n=0:10)时，画出x(n)的频谱X(k) 的幅度；

（2） 将（1）中的x(n)以补零的方式，使x(n)加长到(n:0~100)时，

画出x(n)的频谱X(k) 的幅度；

（3） 取x(n)(n:0~100)时，画出x(n)的频谱X(k) 的幅度。利用FFT

进行谱分析 已知：模拟信号

以t=0.01n(n=0:N-1)进行采样，求N 点DFT 的幅值谱。 请分别画出N=45; N=50;N=55;N=60时的幅值曲线。

数字信号处理实验一

1.(1) L=5;N=20;

60

,7)4(;60,5)3(;40,5)2(;20,5)1()]

(~[)(~

,2,1,01)1(,01,1)(~=========±±=??

?-+≤≤+-+≤≤=N L N L N L N L n x DFS k X m N m n L mN L mN n mN n x )

52.0cos()48.0cos()(n n n x ππ+=)

8cos(5)4sin(2)(t t t x ππ+=

n=1:N;

xn=[ones(1,L),zeros(1,N-L)];

Xk=dfs(xn,N);

magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]);

k=[-N/2:N/2];

figure(1)

subplot(2,1,1);

stem(n,xn);xlabel('n');ylabel('xtide(n)'); title('DFS of SQ.wave:L=5,N=20');

subplot(2,1,2);

stem(k,magXk);

axis([-N/2,N/2,0,16]);

xlabel('k');ylabel('Xtide(k)');

(2)

L=5;N=40;

n=1:N;

xn=[ones(1,L),zeros(1,N-L)];

Xk=dfs(xn,N);

magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]);

k=[-N/2:N/2];

figure(2)

subplot(2,1,1);

stem(n,xn);xlabel('n');ylabel('xtide(n)'); title('DFS of SQ.wave:L=5,N=40');

subplot(2,1,2);

stem(k,magXk);

axis([-N/2,N/2,0,16]);

xlabel('k');ylabel('Xtide(k)');

(3)

L=5;N=60;

n=1:N;

xn=[ones(1,L),zeros(1,N-L)];

Xk=dfs(xn,N);

magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]);

k=[-N/2:N/2];

figure(3)

subplot(2,1,1);

stem(n,xn);xlabel('n');ylabel('xtide(n)'); title('DFS of SQ.wave:L=5,N=60');

subplot(2,1,2);

stem(k,magXk);

axis([-N/2,N/2,0,16]);

xlabel('k');ylabel('Xtide(k)');

(4)

L=7;N=60;

n=1:N;

xn=[ones(1,L),zeros(1,N-L)];

Xk=dfs(xn,N);

magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]);

k=[-N/2:N/2];

figure(4)

subplot(2,1,1);

stem(n,xn);xlabel('n');ylabel('xtide(n)'); title('DFS of SQ.wave:L=7,N=60');

subplot(2,1,2);

stem(k,magXk);

axis([-N/2,N/2,0,16]);

xlabel('k');ylabel('Xtide(k)');

2. (1)

M=10;N=10;n=1:M;

xn=cos(0.48*pi*n)+cos(0.52*pi*n);

n1=[0:1:N-1];y1=[xn(1:1:M),zeros(1,N-M)]; figure(1)

subplot(2,1,1);stem(n1,y1);xlabel('n'); title('signal x(n),0<=n<=10');

axis([0,N,-2.5,2.5]);

Y1=fft(y1);magY1=abs(Y1(1:1:N/2+1));

k1=0:1:N/2;w1=2*pi/N*k1;

subplot(2,1,2);

title('Samples of DTFT Magnitude');

stem(w1/pi,magY1); axis([0,1,0,10]);

xlabel('frequency in pi units');

(2)

M=10;N=100;n=1:M;

xn=cos(0.48*pi*n)+cos(0.52*pi*n);

n1=[0:1:N-1];y1=[xn(1:1:M),zeros(1,N-M)]; figure(2)

subplot(2,1,1);stem(n1,y1);xlabel('n'); title('signal x(n),0<=n<=10');

axis([0,N,-2.5,2.5]);

Y1=fft(y1);magY1=abs(Y1(1:1:N/2+1));

k1=0:1:N/2;w1=2*pi/N*k1;

subplot(2,1,2);

title('Samples of DTFT Magnitude');

stem(w1/pi,magY1); axis([0,1,0,10]);

xlabel('frequency in pi units');

(3)

M=100;N=100;n=1:M;

xn=cos(0.48*pi*n)+cos(0.52*pi*n);

n1=[0:1:N-1];y1=[xn(1:1:M),zeros(1,N-M)]; figure(3)

subplot(2,1,1);stem(n1,y1);xlabel('n'); title('signal x(n),0<=n<=100');

axis([0,N,-2.5,2.5]);

Y1=fft(y1);magY1=abs(Y1(1:1:N/2+1));

k1=0:1:N/2;w1=2*pi/N*k1;

subplot(2,1,2);

title('Samples of DTFT Magnitude');

stem(w1/pi,magY1); axis([0,1,0,10]);

xlabel('frequency in pi units');

3.

figure(1)

subplot(2,2,1)

N=45;n=0:N-1;t=0.01*n;

q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N);

plot(q,abs(y))

stem(q,abs(y))

title('FFT N=45')

%

subplot(2,2,2)

N=50;n=0:N-1;t=0.01*n;

q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N);

plot(q,abs(y))

title('FFT N=50')

%

subplot(2,2,3)

N=55;n=0:N-1;t=0.01*n;

q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N);

plot(q,abs(y))

title('FFT N=55')

%

subplot(2,2,4)

N=16;n=0:N-1;t=0.01*n;

q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N);

plot(q,abs(y))

title('FFT N=16')

function[Xk]=dfs(xn,N)

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

k=[0:1:N-1];

WN=exp(-j*2*pi/N);

nk=n'*k;

WNnk=WN.^nk;

Xk=xn*WNnk;

实验二 用双线性变换法设计IIR 数字滤波器 一、 实验目的

1． 熟悉用双线性变换法设计IIR 数字滤波器的原理与方法； 2． 掌握数字滤波器的计算机仿真方法；

3．

通过观察对实际心电图的滤波作用，获得数字滤波器的感性知

识。

二、 实验内容

1．

用双线性变换法设计一个巴特沃斯低通IIR 滤波器，设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB ；在阻

带内[0.3π，π]频率区间上，最小衰减大于1 5dB 。 2． 以0.2π为采样间隔，打印出数字滤波器在频率区间[0，0.2π]上的幅值响应曲线。

3．

用所设计的滤波器对实际的心电图信号采样序列x （n ）（见教科书上254页实验二的第六部分）进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。

三、 实验步骤

1．（参考数字信号处理教材实验二）；

2．（参考数字信号处理教材实验二）；

3．用MATLAB 语言编写仿真滤波程序，完成实验内容1~3； 4．写出完整的实验报告并回答教科书上的思考题。

四、 I IR 滤波器设计的常用函数

1. Butter

功能：Butterwoeth （巴特沃斯）数字滤波器设计。 格式：[b ，a]=butter （n ，Wn ）

可设计出截止频率为Wn 的n 阶低通Butterwoeth 滤波器，其滤波器为

n

n

z

n a z

a z

n b z b b z A z B z H ----++++++++=

=

)1()2(1)1()2()1()

()()(11

截止频率是滤波器幅度下降至2

/

1处的频率.

2．Buttord

IIR(巴特沃斯)滤波器阶的选择

格式 [n,Wn]=buttord(Wp,Ws,Rp,Rs)----------数字域

[n,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)-----模拟域

说明：buttord 可在给定滤波器性能的情况下，选择模拟或数字Butterword 滤波器的最小阶，其中Wp 和Ws 分别是通带和阻带的截止频率，其值

1)(0≤≤s

p orW

W ，

当其值为1时，表示0.5S f ,R p ,R S 分别是通带和阻带区的波纹系数。

例如：设计一低通滤波器，通带范围0~100Hz ，通带波纹小于3dB ，阻带为-30dB ，

并利用最小的阶来实现，其程序如下

Wp=100/500； Ws=200/500； [n ，Wn]=buttord （Wp,Ws,3，30）； [b ，a]=butter （n ，Wn ）；

2-1、 T=pi/2;Fs=1/T ;

wpz=0.2;wsz=0.3;

wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=1;rs=15; [N,wc]=buttord(wp,ws,rp,rs,'s') [B,A]=butter(N,wc,'s')

k=0:511;

fk=0:0.2*pi:pi/2;

wk=2*pi*fk;

Hk=freqs(B,A,wk);

subplot(1,1,1);

plot(fk,20*log10(abs(Hk)));grid on

xlabel('频率（Hz）');ylabel('幅度（dB）')

2-2

[b,a]=butterw(0.2,0.3,1,15); figure(1)

freqz(b,a,[0:0.02*pi:0.2*pi])

[b,a]=butterw(0.2,0.3,1,15);

xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38 ,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4 ,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

N=length(xn);

figure(2)

subplot(2,1,1);plot(xn)

title('before')

yn=filter(b,a,xn);

subplot(2,1,2);plot(yn)

title('after')

准备函数：

Butterw

function [b,a]=butterw(wp,ws,ap,as) T=0.2;

Wp=(2*tan(wp/2))/T;

Ws=(2*tan(ws/2))/T;

t1=0.1*ap;

t2=0.1*as;

k=((10^t1-1)/(10^t2-1))^0.5;

r=Ws/Wp;

n=-log10(k)/log10(r);

n=floor(n)+1;

Wn=Wp*(10^0.1*ap-1)^(-1/(2*n)); [Z,P,K]=butter(n,Wn,'s');

lp=length(P);

z=-ones(1,lp);

z=z';

p=1/(1-P);

k=K/((2/T)^lp);

[b,a]=zp2tf(z,p,k);

实验三 用窗函数法设计FIR 数字滤波器

1. 实验目的及原理方法参考（301~302）页

2. 实验内容

● 用MATLAB 产生各种窗函数 ● 利用窗函数设计FIR 滤波器 参考程序：

（1）用MATLAB 实现各种窗函数

m=200;

a=boxcar(m); 矩形窗 b=hanning(m); 汉宁窗 c=hamming(m); 海明窗

d=blackman(m); 布莱克曼窗 e=kaiser(m,7.865); 凯泽窗 m=1:200;

plot(m,a,'r*',m,b,'g+',m,c,'y*',m,d,'b.')

（2） 利用窗函数设计FIR 滤波器

设计具有下列指标ωp =0.25π,Rp=0.25dB ，ωs =0.3π,Rp=50dB 的低通数字滤波器。 要求：

● 选择合适的窗函数; ● 画出滤波器的频率特性

几个主要的MATLAB 函数

（1）[r ,p,k]=residuez(b,a) 将有理分式变换成部分分式表示 （2）[b,a]=residuez(r ,p,k) 将部分分式换成有理分式变表示

式中各参数意义如下：

n

n m m z

a z

a a z

b z

b b z H ----++++++=

1

10110)()

(1

1)1()1()(1)()1(1)

1(n m z n m k k z

n p n r z p r +---+-++-+-=

（3）freqz(b,a) 画连续系统的频率特性曲线

（4）impz(b,a，101)画系统的脉冲响应

三、实验程序

1}

m=200;

a=boxcar(m);%矩形窗

b=hanning(m);%汉宁窗

c=hamming(m);%海明窗

d=blackman(m);%布莱克曼窗

e=kaiser(m,7.865);%凯泽窗

m=1:200;

plot(m,a,'r*',m,b,'g+',m,c,'y*',m,d,'b.')

2)

Wp=0.25*pi;

Ws=0.3*pi;

tr_width=Ws-Wp;%求过渡带宽

N=ceil(6.6*pi/tr_width)+1%求窗的宽度

n=0:1:N-1;

Wc=(Ws+Wp)/2;%理想低通滤波器的截止频率

hd=ideal_lp(Wc,N);%理想低通滤波器的单位脉冲响应w_ham=(hamming(N))';

h=hd.*w_ham;

[db,mag,pha,grd,w]=freqz_m(h,[1]); delta_w=2*pi/1000;

Ap=-(min(db(1:1:Wp/delta_w+1)))

As=-round(max(db(Ws/delta_w+1:1:501))) subplot (221)

stem (n,hd)

title ('理想单位脉冲响应hd（n）') subplot (222)

stem (n,w_ham)

title ('海明窗w (n)')

subplot (223)

stem (n,h)

title ('实际单位脉冲响应hd (n)')

subplot (224)

plot (w/pi,db)

title ('幅度响应（dB）')

axis ([0,1,-100,10])

准备函数

Freqz

function [db,mag,pha,grd,w]=freqz_m(b,a); [H,w]=freqz(b,a,1000,'whole');

H=(H(1:501))';w=(w(1:501))';

mag=abs(H);

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

pha=angle(H);

grd=grpdelay(b,a,w);

ideal

function hd=ideal_lp(wc,M);

alpha=(M-1)/2;

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

m=n-alpha+eps;

hd=sin(wc*m)./(pi*m);