MatLab仿真试验报告
2008211208班
08211106号
史永涛
班内序号：01
实验一：数字信号的 FFT 分析
一． 实验目的
通过本次实验，应该掌握：
(1) 用傅立叶变换进行信号分析时基本参数的选择。
(2) 经过离散时间傅立叶变换（DTFT）和有限长度离散傅立叶变换（DFT） 后信号频谱上的区别，前者 DTFT 时间域是离散信号，频率域还是连续的，而 DFT 在两个中都是离散的。
(3) 离散傅立叶变换的基本原理、特性，以及经典的快速算法（基2时间抽选法），体会快速算法的效率。
(4) 获得一个高密度频谱和高分辨率频谱的概念和方法，建立频率分辨率和时间分辨率的概念，为将来进一步进行时频分析（例如小波）的学习和研究打下基础。
(5) 建立 DFT 从整体上可看成是由窄带相邻滤波器组成的滤波器组的概念，此概念的一个典型应用是数字音频压缩中的分析滤波器，例如 DVD AC3 和MPEG Audio。
对于第一个小题目
二. 实验任务
1、实验内容及要求
(1) 离散信号的频谱分析：
设信号
此信号的0.3pi 和 0.302pi两根谱线相距很近，谱线 0.45pi 的幅度很小，请选择合适的序列长度 N 和窗函数，用 DFT 分析其频谱，要求得到清楚的三根谱线
三．实验分析
(1)可先做图，得到时域序列的图像
(2)要得到三根谱线，则先用FFT得到时域信号x（n）的频谱X(N),
(3)对时域做图时，应使用离散傅立叶变换，才可得到相互离散的谱线
(4) 要得到题目要求的三根谱线，应选择合适的窗函数的宽度与幅度，否则无法得到0.45pi处的谱线，同时进行变换时，要选择合理的fft点数，否则无法将0.3pi 和 0.302pi两根谱线相分离。
(5)我没有直接得角频率，而是通过式子w=2*pi/k*[0:1:k-1]，得到角频率，然后根据ω = ΩT，（T=1/Fs），可得到Ω和f,从而得到最后的结果。
四．实验代码及波形图
k=1000;
n=[1:1:k];
x=0.001*cos(0.45*n*pi )+ sin(0.3*n*pi )- cos(0.302*n*pi-pi/4 );
subplot(2,1,1);
stem(n,x,'.');
title('时域序列');
xlabel('n');
ylabel('x(n)');
xk=fft(x,k);
w=2*pi/k*[0:1:k-1];%数字角频率
subplot(2,1,2);
stem(w/pi,abs(xk));
axis([0 0.5 0 2]);%关键代码设置窗函数的宽度与限幅
title('1000点dft');
xlabel('数字频率');
ylabel('|xk(k)|');
所得实验图像为 ：
如图，得到了清晰的0.3pi, 0.302pi,0.45pi处的谱线 ，完成了实验要求。
对于第二个小题目：
二 .实验任务
DTMF 信号频谱分析：
用计算机声卡采用一段通信系统中电话双音多频（DTMF）拨号数字 0～9的数据，采用快速傅立叶变换（FFT）分析这10个号码DTMF拨号时的频谱。
三．实验分析
（1）对于题目要求的用计算机声卡采集0-9的DTMF，可使用第二题中的到的自己手机号码的DTMF
（2）画出所得脉冲DTMF信号的时域波形
（3）用离散傅立叶变换，得到信号的频谱
（4）同样应注意选择合适的FFT点数，以及窗函数的宽度和幅度。
四．实验代码和波形图
[x,Fs,bits]=wavread('syt.wav');
N=1:1:Fs;
subplot(2,1,1);
plot(x);
title('时域波形');
xlabel('n');
ylabel('x(n)');
n=Fs*2
xk=fft(x,n);
k=0:1:n-1;
w=2*pi/n*k;
subplot(2,1,2);
stem(w/pi,abs(xk),'.');
title('音频信号频谱');
xlabel('数字频率');
axis([0 2 0 1000])
所得实验图像为
如图，得到了一个总的频谱图，完成了实验的任务要求。
五．实验中遇到的问题及解决方法
（1），在得到三根清晰的频谱时，不知道如何选择合适的FFT取样点数，只好不断进行取点实验，花费了不少时间，最后还是通过询问同学以及上网查找资料，才得到了最佳取样点,N=1000。
（2）在分析这10个号码DTMF拨号时的频谱，不会使用wavread函数，老是弄错等号左边的格式，后经查询资料才学会使用wavread函数。
六．心得体会
在实际的数字系统中，DFT是一种得到了广泛应用的、重要的信号处理手段。其之所以重要，不仅仅得益于DFT能够更好地适应于数字系统的两个基本特征并拥有严格的定义和明确的物理含义，同时还在于其具备了快速计算方法，为FFT（fast fourier transform），FFT极大地降低了DFT的运算量。因此我们应当了解DFT的原理，掌握FFT的使用方法，同时应会利用matlab软件求解较为简单的FFT，从而分析一些信号的频谱掌握傅立叶变化的原理特性，为以后的学习打下基础。
本次试验是在离散信号时域波形的基础上，对离散信号做傅里叶变换而得到频域波形进行分析。在进行变换时要合理的选择fft点数，以及窗函数的宽度与幅度。通过本次试验，我基本掌握了用傅立叶变换进行信号分析时基本参数的选择，离散傅立叶变换的基本原理、特性，以及经典的快速算法，为以后dsp及matlab的学习打下了基础。这个实验难度虽然比后两个实验略小，但由于之前没有实际使用过matlab仿真dsp，因此需要先从头学起，掌握一些常用的函数和操作，从而也为之后的两个实验打下了基础。
实验二 DTMF 信号的编码
一、实验目的
(1)复习和巩固 IIR 数字滤波器的基本概念；
(2)掌握 IIR 数字滤波器的设计方法；
(3)掌握 IIR 数字滤波器的实现结构；
(4)能够由滤波器的实现结构分析滤波器的性能（字长效应）；
(5)了解通信系统电话 DTMF 拨号的基本原理和 IIR 滤波器实现方法。
对于第一个小题目
二． 实验任务
1）把您的联系电话号码 通过DTMF 编码生成为一个 .wav 文件。
（a）根据 ITU Q.23 建议，DTMF 信号的技术指标是：传送/接收率为每秒 10 个号码，或每个号码 100ms。
（b）每个号码传送过程中，信号存在时间至少 45ms，且不多于 55ms，100ms 的其余时间是静音。
（c）在每个频率点上允许有不超过 ±1.5% 的频率误差。任何超过给定频率 ±3.5% 的信号，均被认为是无效的，拒绝接收。
（其中关键是不同频率的正弦波的产生。可以使用查表方式模拟产生两个不同频率的正弦波。正弦表的制定要保证合成信号的频率误差在±1.5%以内，同时使取样点数尽量少）
三． 实验分析
（1）本实验用 IIR 数字滤波器产生双音多频 DTMF 拨号信号的产生和恢复，演示 IIR 数字滤波器在通信系统中的应用。DTMF 信号是将拨号盘上的 0~F 共16 个数字，用音频范围的 8 个频率来表示的一种编码方式。8 个频率分为高频群和低频群两组，分别作为列频和行频。
每个字符的信号由来自列频和行频的两个频率的正弦信号叠加而成。频率组合方式如下图所示。
频率 1209Hz 1336Hz 1477Hz 1633Hz
697Hz 1 2 3 A
770Hz 4 5 6 B
852Hz 7 8 9 C
941Hz * 0 # D
则以此为基础，可得到各数字的DTMF。
（2） 通过设置占空比，使传送号码时，信号存在的时间满足题目的要求。
（3）可以使用查表方式模拟产生两个不同频率的正弦波。正弦表已满足对于频率误差的要求
四．实验代码及波形图
N=800;
total_x=[];
tm=[49,50,51,65;52,53,54,66;55,56,57,67;42,48,35,68];
n=1:N;%用于设置占空比
fs=8000;%用于频率变换
lf=[697 770 852 941];
hf=[1209 1336 1477];
x1=sin(2*pi*lf(1)*n/fs)+sin(2*pi*hf(1)*n/fs);%1
x2=sin(2*pi*lf(1)*n/fs)+sin(2*pi*hf(3)*n/fs);%3
x3=sin(2*pi*lf(2)*n/fs)+sin(2*pi*hf(2)*n/fs);%5
x4=sin(2*pi*lf(3)*n/fs)+sin(2*pi*hf(2)*n/fs);%8
x5=sin(2*pi*lf(1)*n/fs)+sin(2*pi*hf(1)*n/fs);%1
x6=sin(2*pi*lf(2)*n/fs)+sin(2*pi*hf(2)*n/fs);%5
x7=sin(2*pi*lf(1)*n/fs)+sin(2*pi*hf(2)*n/fs);%2
x8=sin(2*pi*lf(2)*n/fs)+sin(2*pi*hf(3)*n/fs);%6
x9=sin(2*pi*lf(1)*n/fs)+sin(2*pi*hf(2)*n/fs);%2
xa=sin(2*pi*lf(3)*n/fs)+sin(2*pi*hf(1)*n/fs);%7
xb=sin(2*pi*lf(2)*n/fs)+sin(2*pi*hf(1)*n/fs);%4
total_x=[x1,x2,x3,x4,x5,x6,x7,x8,x9,xa,xb];
x11=[x1,zeros(1,800)];
x22=[x2,zeros(1,800)];
x33=[x3,zeros(1,800)];
x44=[x4,zeros(1,800)];
x55=[x5,zeros(1,800)];
x66=[x6,zeros(1,800)];
x77=[x7,zeros(1,800)];
x88=[x8,zeros(1,800)];
x99=[x9,zeros(1,800)];
xaa=[xa,zeros(1,800)];
xbb=[xb,zeros(1,800)];
x=[x11,x22,x33,x44,x55,x66,x77,x88,x99,xaa,xbb];
x=x/max(abs(x));
plot(x);
sound(x);
wavwrite(x,fs,'syt.wav')
所得图像为：
如图，得到了11个数字DTMF编码的时域波形，完成了题目的要求。
对于第二个小题目
二．实验任务
对所生成的DTMF文件进行解码。
（1）DTMF 信号解码可以采用 FFT 计算 N 点频率处的频谱值，然后估计出所拨号码。但 FFT计算了许多不需要的值，计算量太大，而且为保证频率分辨率，FFT的点数较大，不利于实时实现。因此，FFT 不适合于 DTMF 信号解码的应用。
（2）由于只需要知道 8 个特定点的频谱值，因此采用一种称为 Goertzel 算法的 IIR 滤波器可以有效地提高计算效率。其传递函数为：
三．实验分析
（1）解码的过程的核心是使用Goertzel 算法，由于只有8种可能的音频要检测，所以需要8个算式给出的滤波器，每个滤波器调谐到这8个频率之上，在相应的频率上的频谱值最大，可设置一个门限来判断，大于门限才予以正确，判断出在频率表中的行和列，在找出对应的ASCII码。根据人们的经验，对DTMF信号最205点的FFT最好，而且各个频率对应的K值都已经算好了（如K=18对应696hz；K=20对应770hz；K=22对应852HZ）；则根据频谱图上的在K为多少时频率最大，则对应出相应的频率，在对应出相应的数字。
（2）使用Goertzel 算法画出每个声音信号对应的8个频率，得到其中最大的两个频率，即一个为行频，一个为列频。
(3)将这两个频率与原频率方式表进行对比，确定这两个频率与表中哪两个频率最为接近，从而得到由这两个频率确定的一个数字。
(4)编码过程中，可以使用一个pause函数，使解码得到的频谱和数字逐一显示，使过程更为清晰。
四．实验代码及波形图
k = [18 20 22 24 31 34 38 42]; % 要求的DFT样本序号
N=205;
disp(['接收端检测到的号码为'])
for nn=1:11
m=800*(nn-1);
X=goertzel(total_x(m+1:m+N),k+1);% 用Goertzel算法计算八点DFT样本
val = abs(X); % 列出八点DFT向量
stem(k,val,'.');
grid;
xlabel('k');ylabel('|x(k)|')
set(gcf,'color','w')
shg;
pause;
limit = 80; %
%
for s=5:8;
if val(s) > limit, break, end % 查找列号
end
for r=1:4;
if val(r) > limit, break, end % 查找行号
end
disp([setstr(tm(r,s-4))]) % 显示接收到的字符
end
所得到的图像为
此为接收到第一个声音信号后解码得到的频谱，如图，分析其最接近的行频和列频，与原表进行对比，可得到这两个频率对应的数字即为“1”。
此时程序显示为：
接收端显示的号码为：1
此为接收到第二个声音信号后解码得到的频谱，如图，分析其最接近的行频和列频，与原表进行对比，可得到这两个频率对应的数字即为‘3’。
此时程序显示为：
接收端显示的号码为：3
其他数字的接收显示与上两个数字原理相同，就不在此一一列举了，
最终显示为：
五．实验中遇到的问题及解决方法
（1）对于第一个小题目，我采用的是穷举法，得到自己电话号码中每一个数字对应的两个频率，可以得到DTMF编码，但没有一个方便的人机交互界面，且只适用于自己的电话号码。
（2）在解码时，由于在程序中放错了pause的位置，使得程序的输出与图像的显示顺序出了点小问题，后经查书及时得到了改正。
六．心得体会
首先，本次实验通过查表自己使用正弦信号产生单音多频的声音信号，从而得到一串号码的声音信息。通过实验了解了单音多频的原理及应用。不过，值得注意的是，实验中要考虑对于音频信号占空比及频率变换的选择。
然后，在解码过程中，要了解Goertzel 算法，由于只有8种可能的音频要检测，所以需要8个算式给出的滤波器，每个滤波器调谐到这8个频率之上，在相应的频率上的频谱值最大，可设置一个门限来判断，大于门限才予以正确，这样就得到没个信号对应的数字，最终解出一个完整的电话号码。
这次的实验题目明显比上一个难度有所增加，但由于我们对于matlab软件有了些许的了解，所以做起来也不是特别困难，尤其是这个实验题目，是用我们自己电话号码进行编码和解码，更是增加了我们的兴趣，通过这次实验，使我们对于DTMF编码和解码有了些许的了解，同时进一步加深了我们对于matlab软件的使用。
实验三：FIR 数字滤波器的设计和实现
一． 实验目的
（1）掌握FIR 数字滤波器窗口设计法的原理和设计步骤；
（2）了解Gibbs 效应发生的原因和影响；
（3）了解不同类型的窗函数对滤波效果的影响，以及窗函数和长度 N 的选择。
二． 实验任务
（1）录制自己的一段声音，长度为 45秒，取样频率 32kHz，然后叠加一个高斯白噪声，使得信噪比为 20dB。请采用窗口法设计一个 FIR 通滤波器，滤除噪声提高质量。
（2）滤波器指标参考：通带边缘频率为 4kHz，阻带边缘频率为4.5kHz，阻带衰减大于 50dB；
（3）Matlab 函数 y = awgn(x,snr,'measured') ，首先测量输入信号 x 的功率，然后对其叠加高斯白噪声；
三． 实验分析
（1） 滤波器通带边缘为4KHZ，组带边缘频率为4.5Khz，可得实为低通滤波器
（2） 滤波器带通范围为0～4Khz，则要通过滤波器的声音信号必为人的声音。
（3） 阻带衰减大于50db，则查书得滤除噪声，且满足此要求，使用汉明窗较为合适
（4） 画出原始信号，加噪声后以及滤波后信号的频谱，进行对比，观察滤波器的性能
（5） .最终用wavewrite函数，得到加噪声和滤波后的声音信号从另一个方面进行对比
（6） 画出滤波器的幅频相应和相频相应，观察滤波器的性能。
四． 实验代码和波形图
%----------------------------------------------1，得到原始信号的频谱
[y1,fs,nbits]=wavread('mal.wav');
fs=38000;
k=1:4096;
Yk1=fft(y1,4096);
fprintf('按回车输出原始信号频谱\n');
pause;
figure(1)
subplot(2,1,1);
plot(y1);
title('原始文件时域');
subplot(2,1,2);
plot(32/4096*k,abs(Yk1))
axis([0,6,0,15]);
xlabel('f/KHZ')
title('原始文件频域');
%---------------------------------------------2，叠加噪声后的频谱
y2= AWGN(y1,20,'measured','db');
fprintf('按回车输出加噪声后频谱\n');
pause;
Yk2=fft(y2,4096);
figure(2)
subplot(2,1,1);
plot(y2);
title('加噪声后的时域');
subplot(2,1,2);
plot(32/4096*k,abs(Yk2))
xlabel('f/KHZ');
axis([0,6,0,15]);
title('加噪声后的频域');
%------------------------------------------3，使用filter函数设计滤波器
fp=4000;
fr=4500;
wp=2*pi*fp/fs;
wr=2*pi*fr/fs;
tr_width=wr-wp;
N=ceil(6.6*pi/tr_width)+1;
n=0:1:N;
wc=(wr+wp)/2;
alpha=(N-1)/2;
n=0:1:N-1;
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);
b=fir1(N,wc/(4*pi/3.65));
w_ham=(hamming(N))';
h=hd.*w_ham;
[H,w]=freqz(h,[1],1000,'whole');
H=(H(1:1:501))';
w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
delta_w=2*pi/1000;
x=filter(b,1,y2);
%-------------------------------------------4，使用得到的滤波器进行滤波
fprintf('按回车输出滤波后频谱\n');
pause;
figure(3)
subplot(2,1,1);
plot(x);
title('滤除噪声后的输出');
Xk=fft(x,4096);
subplot(2,1,2);
plot(32/4096*k,abs(Xk))
axis([0,6,0,15]);
xlabel('f/KHZ')
fprintf('按回车输出滤波器频谱\n');
pause;
%--------------------------------------------5，画出滤波器的频谱
wn=(ws+wp)/2;
t=0:1/fs:(size(x2)-1)/fs; %选择窗函数，并归一化截止频率，
figure(4) % 得滤波器相应图像
所得到的图像为：
（1） 原始信号时域、频域波形
（2） 加噪声的时域、频域波形为
如图，频域图像中可明显看出加入了噪声信号。
（3） 滤除噪声后的时域、频域图像
如图，可以看出，在4Khz以后，滤波器作用明显，噪声被滤除掉了，符合题目要求。
（4） 最终得到的滤波器的幅频和相频特性图像
、
五． 实验中遇到的问题及解决方法
（1）开始时，不会选择窗函数，由于理论课程还没涉及滤波器的知识，所以必须自己自学滤波器的相关知识，从而进行设计。
（2）得到滤波器后，本来选择通过卷积得到滤波后的信号，可是根本不能通过编译，后来不得不重新选择使用filter函数，才得到了滤波后的输出。
（3）对于输出的滤波器的频谱图像不是很理解，后经查书得知其为归一化的幅度，相位谱，故横坐标都限制在0～1内。
六． 心得体会
这个题目综合性较强，主要是针对滤波器的设计与使用，让我们了解了滤波器的工作原理和使用方法，了解了各种窗函数的选择与使用为今后的理论课学习打下了良好的基础。同时涉及到了许多新的matlab中的函数，
让我们增强了对于matlab的使用，有助于我们在今后的学习过程中更加得心应手的使用matlab这一强大的辅助工具。
四．实验总结
短短的三次matlab仿真实验已经结束了，我感觉自己从中收获良多，不仅仅是对于知识的理解，对于软件的掌握，更多是一种学习方法和态度。
通过这几次课程，我们学会了用matlab显示出FFT的频谱波形，学会了IIR和FIR滤波器的设计，学会了DTMF编码和解码，通过这些设计，一方面加深了我们对于课本知识的理解，有助我们今后对于理论知识的进一步深，另一方面加深了我们对于matlab的理解，提高了我们的编程能力。
在进行仿真的过程中，我的确遇到了不少的问题，首先就是理论知识跟不上，由于课本还没学到这，我们不得不自学来进行最初步的理解；然后就是对于matlab这门语言掌握的太差，刚开始时，连一些基本的函数，都缺乏了解，后来经过几次询问老师，和课下进行查书，才慢慢摸着了门道，最终得到了正确的结果。
通过这次实验，我感觉自己最大的收获就是掌握了一种新的学习方法：用软件进行仿真，同时于理论知识相结合，既能增加自己对于理论知识的理解，又能提高自己的动手能力，一举多得。
最后，感谢老师对于我们的教导，让我们收获良多。