matlab仿真实验报告

MATLAB仿 真 实 验 报 告

院 系：电子工程学院

姓名：王 力

班 级： 2013211207

学 号： 2013211006

实验一：数字信号的 FFT 分析

实验内容及要求

离散信号的频谱分析：

设信号 x( n) 0.001*cos( 0.45n ) sin( 0.3n ) cos(0.302n )

4

此信号的 0.3pi 和 0.302pi 两根谱线相距很近，谱线 0.45pi 的幅度很小，请选择合适的序列长度 N 和窗函数，用 DFT 分析其频谱， 要求得到清楚的三根谱线 。

DTMF 信号频谱分析 P218 -225 4.9.3 双音） 用计算机声卡采用一段通信系统中

电话双音多频（ DTMF）拨号数字 0～9 的数据，采用快速傅立叶变换（ FFT）分析这 10 个号码 DTMF 拨号时的频谱。

实验分析

要得到清晰的三根谱线，用 matlab 内置函数 fft 对时域信号进行快速傅里叶变换，需要选好

变换点数 N ，以避免出现频谱模糊现象。程序中选择 N=1000

由于谱线 0.45pi 的幅度很小，在作图时需要对坐标比例进行控制。使用 axis 函数实现。

1. 代码及注释

频谱分析 ：

N = 1000; % Length of DFT

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

xn = 0.001*cos(0.45*n*pi)+sin(0.3*n*pi)-cos(0.302*n*pi-pi/4);

Xk = fft(xn,N);

xn1=0.001*cos(0.45*n*pi);

Xk1 = fft(xn1,N);

xn2=sin(0.3*n*pi);

Xk2 = fft(xn2,N);

xn3=-cos(0.302*n*pi-pi/4);

Xk3 = fft(xn3,N);

k=[0:1:N/2];

w=2*pi*k/N;

subplot(2,1,1);

stem(w/pi,abs(xn(1:1:(N/2+1))),'.');

title('x(n)');xlabel('w/pi');

axis([0.28,0.50,0,4]);

subplot(2,1,2);stem(w/pi,abs(Xk(1:1:(N/2+1))),'r.');

title('DFT 0.001*cos(0.45*n*pi)+sin(0.3*n*pi)-cos(0.302*n*pi-pi/4)'); xlabel('w/pi');axis([0.28,0.50,0,4]);

结果截图：

DTMF 频谱分析：

clear; close all;

column=[1209,1336,1477,1633];

line=[697,770,852,941];

fs=10000;

N=1024;

ts=1/fs;

n=0:N-1;

f=0:fs/N:fs/N*(N-1);

key=zeros(16,N);

key(1,:)=cos(2*pi*column(1)*n*ts)+cos(2*pi*line(1)*n*ts);

key(2,:)=cos(2*pi*column(2)*n*ts)+cos(2*pi*line(1)*n*ts);

key(3,:)=cos(2*pi*column(3)*n*ts)+cos(2*pi*line(1)*n*ts);

key(4,:)=cos(2*pi*column(1)*n*ts)+cos(2*pi*line(2)*n*ts);

key(5,:)=cos(2*pi*column(2)*n*ts)+cos(2*pi*line(2)*n*ts);

key(6,:)=cos(2*pi*column(3)*n*ts)+cos(2*pi*line(2)*n*ts);

key(7,:)=cos(2*pi*column(1)*n*ts)+cos(2*pi*line(3)*n*ts);

key(8,:)=cos(2*pi*column(2)*n*ts)+cos(2*pi*line(3)*n*ts);

key(9,:)=cos(2*pi*column(3)*n*ts)+cos(2*pi*line(3)*n*ts);

key(10,:)=cos(2*pi*column(2)*n*ts)+cos(2*pi*line(4)*n*ts);

figure;

for i=1:10

subplot(4,4,i),plot(f,abs(fft(key(i,:))))

grid;

end

结果截图：

遇到的问题和解决方法

在进行第一题的 DFT 变换时，不知道该选取多大的取样点数

通过计算确定数字信号的最小为 1000 周期才选定了 N=1000 。



N 才能得到清晰的三根谱线。

实验二： DTMF 信号的编码和解码

实验内容及要求

1）把您的联系电话号码 通过 DTMF 编码生成为一个 .wav 文件。

技术指标：

根据 ITU Q.23 建议， DTMF 信号的技术指标是：传送 /接收率为每

秒 10 个号码，或每个号码 100ms。

每个号码传送过程中，信号存在时间至少 45ms，且不多于 55ms，

100ms 的其余时间是静音。

在每个频率点上允许有不超过 ± 1.5% 的频率误差。任何超过给定频率 ± 3.5% 的信号，均被认为是无效的，拒绝接收。

（其中关键是不同频率的正弦波的产生。可以使用查表方式模拟产生两个不同频率的正弦

波。正弦表的制定要保证合成信号的频率误差在± 1.5%以内，同时使取样点数尽量少）

2）对所生成的 DTMF 文件进行解码。

DTMF 信号解码可以采用

FFT 计算 N 点频率处的频谱值，然后估计出所

拨号码。但 FFT 计算了许多不需要的值，计算量太大，而且为保证频率分

辨率， FFT 的点数较大，不利于实时实现。因此，

FFT 不适合于 DTMF

信

号解码的应用。

由于只需要知道 8 个特定点的频谱值， 因此采用一种称为

Goertzel 算法的

IIR 滤波器可以有效地提高计算效率。其传递函数为：

fs=8kHz.

H k ( z)

1 e j 2 k/ N z 1

2cos(2 k / N ) z 1

z

2

1

实验分析

DTMF 信号是将拨号盘上的 0~F 共 16 个数字，用音频范围的

8 个频率来表示的一种编码方

式。

　8 个频率分为高频群和低频群两组，分别作为列频和行频。

每个字符的信号由来自列频和行频的两个频率的正弦信号叠加而成。 频率组合方式如下图所

示。

根据图片即可得到各个数字对应的 DTFM 信号。

通过 zeros 全零矩阵来设置占空比，以达到题目要求。得到信号后，使用 sound 函数来播放

拨号音， writewave 将信号写入声音文件。

在解码时， 使用 Goertzel 算法。滤波器调谐到这 8 个频率之上后， 在相应的频率上的频谱值

最大，通过与标准值的对比找出在 DTMF 图中的行和列，再对应出相应的拨号数字。

查阅资料知， 205 点的 FFT 最佳，而且每个频率对应的 K 值都已经给定（如 K=18 对应 696hz；

K=20 对应 770hz；K=22 对应 852HZ ），则根据频谱图上最大值对应的 K 值，就可以求出相

应的频率，从而比对得出数字。

代码及注释

N=800;fs=8000;% 每个号码 100ms 800 8000

tm=[49,50,51,65;52,53,54,66;55,56,57,67;42,48,35,68];

n=1:N;% 取样点

flow=[697 770 852 941];% 低频

fhigh=[1209 1336 1477 1633];% 高频

x01=sin(2*pi*flow(1)*n/fs)+sin(2*pi*fhigh(1)*n/fs); %1

x02=sin(2*pi*flow(3)*n/fs)+sin(2*pi*fhigh(2)*n/fs); %8 x03=sin(2*pi*flow(2)*n/fs)+sin(2*pi*fhigh(2)*n/fs); %5 x04=sin(2*pi*flow(1)*n/fs)+sin(2*pi*fhigh(1)*n/fs); %1 x05=sin(2*pi*flow(1)*n/fs)+sin(2*pi*fhigh(1)*n/fs); %1 x06=sin(2*pi*flow(1)*n/fs)+sin(2*pi*fhigh(3)*n/fs); %3 x07=sin(2*pi*flow(1)*n/fs)+sin(2*pi*fhigh(3)*n/fs); %3 x08=sin(2*pi*flow(3)*n/fs)+sin(2*pi*fhigh(3)*n/fs); %9 x09=sin(2*pi*flow(1)*n/fs)+sin(2*pi*fhigh(2)*n/fs); %2 x10=sin(2*pi*flow(2)*n/fs)+sin(2*pi*fhigh(2)*n/fs); %5 x11=sin(2*pi*flow(3)*n/fs)+sin(2*pi*fhigh(1)*n/fs); %7 x=[x01,x02,x03,x04,x05,x06,x07,x08,x09,x10,x11];% 组成矩阵

x01_z=[x01,zeros(1,800)];% 补零

x02_z=[x02,zeros(1,800)];

x03_z=[x03,zeros(1,800)];

x04_z=[x04,zeros(1,800)];

x05_z=[x05,zeros(1,800)];

x06_z=[x06,zeros(1,800)];

x07_z=[x07,zeros(1,800)];

x08_z=[x08,zeros(1,800)];

x09_z=[x09,zeros(1,800)];

x10_z=[x10,zeros(1,800)];

x11_z=[x11,zeros(1,800)];

x_z=[x01_z,x02_z,x03_z,x04_z,x05_z,x06_z,x07_z,x08_z,x09_z,x10_z,x11_z];

x_z=x_z/max(abs(x_z));

subplot(2,1,1);

plot(x_z);

sound(x_z);

audiowrite('mjl_num.wav',x_z,fs);% 写入声音文件

k=[18 20 22 24 31 34 38 42];

N=205;

subplot(2,1,2);

xk=fft(x_z);

mxk=abs(xk);

disp([' 解码得到的号码是： '])

for i=1:11

m=800*(i-1);

X=goertzel(x(m+1:m+N),k+1);%goertzel 算法做变换

v=abs(X);% 求模

xk2=v.^2;

stem(k,v,'.');% 画脉冲图

grid;

xlabel('k');

ylabel('x(k)');

set(gcf,'color','w');

shg;

pause;

limit=80;

xk2

for s=5:8;

if v(s)>limit,break,end

end

for r=1:4;

if v(r)>limit,break,end

end

disp([setstr(tm(r,s-4))])

end

结果截图：

解码得到的号码是：遇到的问题和解决方法

在开始的时候并没有设置占空比， 结果导致拨号音过快。 后来通过在信号后补零来进行占空

比设置。在将信号写入声音文件时，上网查得函数 audiowrite 及其用法。

实验三： FIR 数字滤波器的设计和实现

实验内容及要求

录制自己的一段声音， 长度为 10 秒，取样频率 32kHz ，然后叠加一个高斯白噪声，

使得信噪比为 20dB 。请采用窗口法设计一个FIR 带通滤波器，滤除噪声提高质

量。

提示：

b) 滤波器指标参考：通带边缘频率为 4kHz ，阻带边缘频率为 4.5kHz ，阻带衰减大于

50dB；其实是低通滤波器

c) Matlab 函数 y = awgn(x,snr,'measured') ，首先测量输入信号 x 的功率，然后对其

叠加高斯白噪声；

2、实验目的：

通过本次实验，掌握以下知识：

FIR 数字滤波器窗口设计法的原理和设计步骤；

Gibbs 效应发生的原因和影响；

iii. 不同类型的窗函数对滤波效果的影响，以及窗函数和长度 N 的选择。

实验分析

先确定滤波器的设计指标：



滤波器通带边缘为



4KHZ ，阻带边缘频率为



4.5Khz ，阻带衰减

大于 50db。使用汉明窗。

加噪声时使用 awgn 函数。

　使用作图函数分别画出加噪声后的声音信号和去噪声后的声音信号的时域和频域波形。

使用 matlab 库函数 hamming 来进行滤波。再作出滤波后的图形。最后作出滤波器的幅频特性图。

代码及注释

fs=32000;

[x,fs]=audioread('sound.wav');

snr=20; %信噪比为 20dB

x2=awgn(x,snr,'measured','db'); %叠加高斯白噪声 ?

audiowrite('sound_noise.wav',x,fs);

t=0:1/fs:(size(x2)-1)/fs;

wp=8000*pi/32000;

ws=9000*pi/32000;

wdelta=ws-wp;

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

wn=(ws+wp)/2;

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

freqz(b,1,512)

f2=filter(b,1,x2);

title(' 滤波器幅频、相频特性 ');

figure(2)

subplot(2,1,1)

plot(t,x2)

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

subplot(2,1,2)

plot(t,f2)

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

F0=fft(f2,1024);

f=fs*(0:511)/1024;

figure(3)

y2=fft(x2,1024);

subplot(2,1,1)

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

title(' 滤波前频谱 ')

xlabel('Hz');

ylabel(' 幅度 ');

subplot(2,1,2)

F2=plot(f,abs(F0(1:512)));

title(' 滤波后频谱 ')

xlabel('Hz');

ylabel(' 幅度 ');

audiowrite('filterable.wav',f2,fs);

结果截图 ：

遇到的问题和解决方法

在选择窗函数时， 因为理论课还没有学习相关知识， 不知道选择什么类型的窗好， 在阅读课本相关章节后选定 blackman 窗作为滤波器窗函数。在使用 matlab 进行滤波时，不知道如何

操作。在网上查阅相关资料后学习到了使用 matlab 进行滤波的方法。

实验总结

这次实验把 matlab 与数字信号处理的相关知识很好的联系到一起，让我能在实验的过

程中更加深刻的理解数字信号处理的相关内容， 让我觉得受益颇深。 通过这次试验， 使我更

加深刻的理解 FFT、DTMF 以及 FIR 数字滤波器的原理和应用， 我更好的理解了用傅里叶变

换进行信号分析时基本参数的选择、 离散时间傅里叶变换和有限长度离散傅里叶变换后信号

频谱上的区别、 离散傅里叶变化的基本原理、 特性，以及经典的快速算法的效率， 也更好的

掌握了 FIR 数字滤波器的设计方法、实现结构及其性能。

虽然课程很快结束了， 但是在这次实验中我不仅巩固了数字信号处理的内容知识， 也复习了

matlab 仿真软件的运用，为我今后的学习打下了良好的基础。