数字信号实验报告

实验一：快速傅立叶变换的谱分析

一、实验目的：

学会利用matlab中的FFT函数，即进行信号的谱分析。

二、实验题目：

1．已知：

t1=[0:0.001:0.3];

t2=[0.301:0.001:0.6];

t3=[0.601:0.001:0.9];

t4=[0.901:0.001:1.199];

x1=sin(2*pi*100*t1);

x2=sin(2*pi*50*t2);

x3=sin(2*pi*25*t3);

x4=sin(2*pi*10*t4);

信号s1为

t= 0:0.001:1.199;s1= sin(2*pi*100*t)+sin(2*pi*50*t)+sin(2*pi*25*t)+sin(2*pi*10*t);

信号s2为t=[t1,t2,t3,t4];s2= [x1,x2,x3,x4];

信号s3为t=[t1,t2,t3,t4];s3= [x1,x4,x2,x3];

a.编写程序分别画出信号s1,s2,s3的时域波形和幅频图（参考图1，2，3）。

b.观察信号s1,s2,s3的时域波形和频谱图，思考其幅频图的差别及其原因。

c.信号s1,s2,s3的抽样频率fs为多少？对于1200个点的时域离散序列，对其

FFT后仍为长度为1200个点的序列，即周期N=1200，试分析N与离散时间信号频谱的周期fs的对应关系。各频谱图的频率分辨率为多少Hz？

d.若有信号s4的幅频图与s2的幅频图完全相同（如图2所示），问s4的时域

波形和s2是相同的吗，为什么?

实验程序清单：

t1=[0:0.001:0.3];

t2=[0.301:0.001:0.6];

t3=[0.601:0.001:0.9];

t4=[0.901:0.001:1.199];

x1=sin(2*pi*100*t1);

x2=sin(2*pi*50*t2);

x3=sin(2*pi*25*t3);

x4=sin(2*pi*10*t4);

%信号s1的时域图和频域图

t5=0:0.001:1.199;

x5=sin(2*pi*100*t5)+sin(2*pi*50*t5)+sin(2*pi*25*t5)+sin(2*pi*10*t5);

y5=abs(fft(x5));

f=1000*(0:(1/1200):0.5);

figure(1)

subplot(2,1,1),plot(t5,x5);grid on;axis tight;%波形图

title('时域波形');xlabel('t(s)');

subplot(2,1,2),plot(f,y5(1:601));grid on;axis tight;%幅频图

title('幅频图');xlabel('f(Hz)');

%信号s2的时域图和频谱图

t6=[t1,t2,t3,t4];

x6=[x1,x2,x3,x4];

y6=abs(fft(x6));

figure(2)

subplot(2,1,1),plot(t6,x6);grid on;axis tight; %波形图

title('时域波形');xlabel('t(s)');

subplot(2,1,2),plot(f,y6(1:601));grid on;axis tight; %幅频图

title('幅频图');xlabel('f(Hz)');

%信号s3的时域图和频域图

x7=[x1,x4,x2,x3];

y7=abs(fft(x7));

figure(3)

subplot(2,1,1),plot(t6,x7);grid on;axis tight; %波形图

title('时域波形');xlabel('t(s)');

subplot(2,1,2),plot(f,y7(1:601));grid on;axis tight; %幅频图

title('幅频图');xlabel('f(Hz)');

结果分析

a.编写程序分别画出信号s1,s2,s3的时域波形和幅频图（参考下图）。

b.观察信号s1,s2,s3的时域波形和频谱图，思考其幅频图的差别及其原因。

因为信号s1、s2、s3这三个信号中都有频率为100hz，50hz，25hz，10hz的

频率分量，从波形图可以看出s1只是单纯的含有这四种分量，而s2，s3的频谱图相同，除了这四种分量还有其他分量，因为s1中这四种分量一直存在，但s2，s3中只是特定时间段出现。

c. 信号s1,s2,s3的抽样频率fs为多少？对于1200个点的时域离散序列，对其

FFT后仍为长度为1200个点的序列，即周期N=1200，试分析N与离散时间信号频谱的周期fs的对应关系。各频谱图的频率分辨率为多少Hz？

频谱图横坐标所以两相邻点的频率间隔为1000/1200=5/6hz,所以各频谱图的频率分辨率F=5/6HZ.又因为F=1/tp=1/NT=fs/N,因为N=1200,所以s1，s2，s3的抽样频率fs=NF=1200*5/6=1000HZ.

d .若有信号s4的幅频图与s2的幅频图完全相同（如图2所示），问s4的时域波形和s2是相同的吗，为什么?

若有信号s4的幅频图与s2的幅频图完全相同，s4的时域波形和s2可能是不同的，因为幅频特性图只表明了幅度和频率的关系，即对应频率分量的幅值，而其相位特性未知，时域图是幅度和相位图完全一致时才相同。

实验小结

虽然第一次实验看起来比较简单，但是由于平时不怎么接触Matlab，导致很多对Matlab语法以及命令操作不是很熟悉，但是Matlab本省提供的帮助文档很全，所以遇到不会的命令或者函数可以直接调用帮助文档，可以很快的解决问题。其次就是通过实验，反复研究代码，可以加深对快速傅立叶变换的谱分析的理解。

实验二：IIR数字滤波器设计

一、实验目的：

1、熟悉滤波器的基本概念；

2、了解滤波器的分类；

3、熟悉应用matlab设计各种滤波器的方法。

二、实验相关原理：

1.数字滤波器分类

实现方法

有限长冲激响应----FIR滤波器

无限长冲激响应----IIR滤波器

功能

低通（LP）

高通（HP）

带通（BP）

带阻（BS）

2、数字滤波器设计步骤

一、给出所需要的滤波器的技术指标

二、设计一个H(Z)使其逼近所需要的技术指标

三、实现所设计的H(Z)

（ H(Z)为系统的数学模型：传递函数模型、状态方程模型和零极点增益模型

等）

3、IIR数字滤波器设计步骤

一、按一定规则将给出的数字滤波器的技术指标转换为模拟滤波器的技术指标

二、根据转换后的技术指标使用滤波器阶数选择函数，确定最小阶数N和固有频率Wn

三、运用最小阶数N、固有频率Wn产生模拟滤波器原型

四、运用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器函数

1、Buttord(Wp,Ws,Rp,Rs)

butterworth滤波器阶数选择函数,返回符合要求性质的滤波器最小阶数N 以及固有频率Wn

Wp 通带截至频率

Ws 阻带截至频率

Rp 通带衰减（不超过）

Rs 阻带衰减（不小于）

Wp,Ws是归一化频率，范围是[0，1]，对应д弧度归一化处理： Wp（或Ws）/(fs/2)

例：确定数字低通Butterworth滤波器的阶数和固有频率。

要求：wp=500hz,ws=550hz,rp=1db,rs=50db,fs=2000hz

程序清单：

Wp=500;ws=550;rp=1;rs=50;fs=2000 ;

[N,Wn]=buttord(wp/(fs/2),ws/(fs/2),rp,rs)

2、Butter(N,Wn)

设计数字滤波器函数

格式：[B,A]= Butter(N,Wn) [B,A]是滤波器系数矩阵。

如果Wn=[w1,w2]是一个二元向量,则该函数将设计出带通、带阻数字滤波器3、freqz(B，A)

求线性离散时间系统的传输函数的频率响应

格式:

[H,W]= freqz(B，A)

H:频率响应

W:所要计算响应的频率点组成的向量

4、filter()

一维数字滤波函数

格式:

y= filter(B，A,x)

x:待滤波信号

y:滤波以后的信号

x1=wavread('pb8k.wav');%读取语音信号

实验题目：

用双线性变换法设计Butterworth滤波器对加噪后的语音信号进行滤波。绘制频率响应曲线，画出滤波前后的时域图和频谱图。

实验程序清单：

%双线性变换法设计Butterworth滤波器

fs=8000;

x1=wavread('pb8k.wav');

t=(0:length(x1)-1)/8000;

f=fs*(0:1023)/2048;

A1=0.05;A2=0.10;

d=[A1*cos(2*pi*3800*t)+A2*sin(2*pi*3600*t)]';

x2=x1+d;

%用双线性变换法设计一数字巴特沃斯低通器

fp=3200;

fs=3400;

Fs=8000;

Ts=1/Fs;

wp=2*pi*fp*Ts;

ws=2*pi*fs*Ts;

Rp=1;

Rs=15;

wp1=2/Ts*tan(wp/2); %将模拟截止频率进行预畸变

ws1=2/Ts*tan(ws/2);

[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s'); %选择滤波器的最小阶数

[Z,P,K]=buttap(N); %创建butterworth模拟滤波器

%输入N：滤波器的阶数

%输出Z:零点；P：极点；K：增益

[Bap,Aap]=zp2tf(Z,P,K);

[b,a]=lp2lp(Bap,Aap,Wn); %求模拟滤波器的传递函数

[bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换（即求的数字滤波器的系统函数）

[H,W]=freqz(bz,az); %绘制频率响应曲线

figure(1)

plot(W*Fs/(2*pi),abs(H))

grid on;axis tight;

xlabel('频率（Hz）')

ylabel('频率响应')

title('Butterworth')

f1=filter(bz,az,x2);

figure(2)

subplot(2,1,1)

plot(t,x2) %画出滤波前的时域图

grid on;axis tight;

title('滤波前的时域波形');

subplot(2,1,2)

plot(t,f1); %画出滤波后的时域图

grid on;axis tight;

title('滤波后的时域波形');

y3=fft(f1,2048);

figure(3)

y2=fft(x2,2048);

subplot(2,1,1);

plot(f,abs(y2(1:1024))); %画出滤波前的频谱图grid on;axis tight;

title('滤波前的频谱')

xlabel('Hz');

ylabel('幅度');

subplot(2,1,2)

plot(f,abs(y3(1:1024))); %画出滤波后的频谱图grid on;axis tight;

title('滤波后的频谱')

xlabel('Hz');

ylabel('幅度');

实验波形图：

实验小结

对于IIR数字滤波器的设计还真是遇到很多问题，刚开始是找不到音频文件，寻求助教和同学的帮助后，自己慢慢一步一步操作才得到实验结果。主要是刚开始对于IIR数字滤波器原理不是很理解，导致设计滤波器来不够熟练。经过这次实验，可以说从新认识了IIR滤波器原理以及滤波器的设计原理，日后还需多加

练习，熟练掌握IIR数字滤波器的设计。

实验三：综合实验

一、实验目的

对之前的学习内容做到学以致用，通过实验能够设计滤波器并解决具体问题，能够利用数字信号处理语音或图像等实际问题。

二、实验要求及内容

给一段原始的语音信号（可以是自己录制的一段语音），加上一频率为3.8kHz

的高频余弦噪声和频率为3.6kHz的高频正弦噪声（幅度自己可以选择），用窗函

数设计一滤波器（要求最小阻带衰减为50dB）对加噪后的语音信号进行滤波，

画出滤波器的频率响应曲线，画出滤波前后的时域图和频谱图。

需要用到的函数：

fir1 用窗函数设计FIR滤波器的函数

三、实验设计

首先利用窗函数设计出滤波所需要的FIR滤波器，然后对处理前和处理后的语音信号进行时域波形仿真，再分别进行频谱分析，再检查滤波效果。

滤波器设计

1、给出所要设计的FIR模拟滤波器的技术指标

%模拟指标

fp=3200; %通带截止频率

ft=3400; %阻带截止频率

2、计算数字滤波器的数字指标

%数字指标

wp=2*pi*fp*Ts;

ws=2*pi*ft*Ts;

3、根据允许的过渡带宽度及阻带衰减，初步选择窗函数和N值

wdelta=ws-wp; %过渡带宽度

N=ceil(6.6*pi/wdelta); %取整

wn=(0.8+0.85)*pi/2; %截止频率

[bz,az]=fir1(N,wn/pi,hamming(N+1)); %选择窗函数，并归一化截止频率

4、若选用理想低通逼近，理想低通的截止频率

wn=(0.8+0.85)*pi/2; %截止频率

实验程序清单：

fs=8000;

x1=wavread('pb8k.wav');

t=(0:length(x1)-1)/8000;

f=fs*(0:2047)/4096;

A1=0.05;A2=0.10;

d=[A1*cos(2*pi*3600*t)+A2*sin(2*pi*3800*t)]';

x2=x1+d;

%模拟指标

fp=3200; %通带截止频率

ft=3400;%阻带截止频率

Ts=1/fs;

%数字指标

wp=2*pi*fp*Ts;

ws=2*pi*ft*Ts;

wdelta=ws-wp; %过渡带宽度

N=ceil(6.6*pi/wdelta); %取整

wn=(0.8+0.85)*pi/2; %截止频率

[bz,az]=fir1(N,wn/pi,hamming(N+1)); %选择窗函数，并归一化截止频率figure(1)

freqz(bz,az);

grid on;axis tight;

f2=filter(bz,az,x2);

figure(2)

subplot(2,1,1)

plot(t,x2);

grid on;axis tight;

title('滤波前的时域波形');

subplot(2,1,2)

plot(t,f2);

grid on;axis tight;

title('滤波后的时域波形');

y3=fft(f2,4096);

f=fs*(0:2047)/4096;

figure(3)

y2=fft(x2,4096);

subplot(2,1,1);

plot(f,abs(y2(1:2048)));

grid on;axis tight;

title('滤波前的频谱')

xlabel('Hz');

ylabel('幅度');

subplot(2,1,2)

plot(f,abs(y3(1:2048)));

grid on;axis tight;

title('滤波后的频谱')

xlabel('Hz');

ylabel('幅度');

实验波形图：

1、滤波器频谱特性图，幅度响应图和相位响应图，设计出了线性相位。

由上图可知布莱克曼窗比汉明窗得到的滤波器幅度响应衰减更快，但过渡带要比汉明窗宽。

2、语音信号滤波前后的时域波形图。

3、语音信号滤波前后的频域波形图。

由上图可知,用布莱克满窗得到的低通滤波器可以滤除分量为3600hz，3800hz的噪声，满足题目的要求。