北邮数字信号处理MATLAB实验报告.docx
《北邮数字信号处理MATLAB实验报告.docx》由会员分享，可在线阅读，更多相关《北邮数字信号处理MATLAB实验报告.docx（19页珍藏版）》请在冰豆网上搜索。
北邮数字信号处理MATLAB实验报告
数字信号处理软件实验
——MatLab仿真实验报告
学院:
电子工程学院
班级：
2013211202
姓名：
学号：
实验一：
数字信号的FFT分析
1、实验内容及要求
（1）离散信号的频谱分析：
设信号
此信号的0.3pi和0.302pi两根谱线相距很近，谱线0.45pi的幅度很小，请选择合适的序列长度N和窗函数，用DFT分析其频谱，要求得到清楚的三根谱线。
（2）DTMF信号频谱分析
用计算机声卡采用一段通信系统中电话双音多频（DTMF）拨号数字0～9的数据，采用快速傅立叶变换（FFT）分析这10个号码DTMF拨号时的频谱。
2、实验目的
通过本次实验，应该掌握：
（a）用傅立叶变换进行信号分析时基本参数的选择。
（b）经过离散时间傅立叶变换（DTFT）和有限长度离散傅立叶变换（DFT）后信号频谱上的区别，前者DTFT时间域是离散信号，频率域还是连续的，而DFT在两个域中都是离散的。
（c）离散傅立叶变换的基本原理、特性，以及经典的快速算法（基2时间抽选法），体会快速算法的效率。
（d）获得一个高密度频谱和高分辨率频谱的概念和方法，建立频率分辨率和时间分辨率的概念，为将来进一步进行时频分析（例如小波）的学习和研究打下基础。
（e）建立DFT从整体上可看成是由窄带相邻滤波器组成的滤波器组的概念，此概念的一个典型应用是数字音频压缩中的分析滤波器，例如DVDAC3和MPEGAudio。
3.设计思路及实验步骤
1）离散信号的频谱分析：
该信号中要求能够清楚的观察到三根谱线。
由于频率0.3pi和0.302pi间隔非常小，要清楚的显示，必须采取足够大小的N，使得分辨率足够好，至少到0.001单位级，而频率0.45pi的幅度很小，要清楚的观察到它的谱线，必须采取幅度够大的窗函数，使得它的频谱幅度变大一些。
同时还要注意频谱泄漏的问题，三个正弦函数的周期（2pi/w）分别为20,40,1000，所以为了避免产生频谱泄漏（k=w/w0为整数），采样点数N必须为1000的整数倍。
2）DTMF信号频谱分析
双音多频信号中4*4拨号盘中的每一个按键均对应两个频率，一个高频，一个低频，每个数字由两个不同频率的正弦波组成，高频为：
1209HZ,1336HZ,1633HZ,1477HZ,低频为697HZ,770hz，852hz,941hz。
分别用两个数组来装载高频频率和低频频率，然后再分别产生两路频率混叠的信号。
再利用fft频谱分析做出频谱图。
4.实验代码及实验结果
1）离散信号的频谱分析代码
clf;
closeall;%关闭所有图形窗口
N=1000;%DFT点数
n=[1:
1:
N];
x=0.001*cos（0.45*n*pi）+sin（0.3*n*pi）-cos（0.302*n*pi-pi/4）;
y=fft（x,N）;
mag=abs（y）;%对FFT结果求模
w=2*pi/N*[0:
1:
N-1];%数字角频率w
subplot（2,1,1）;%将图形窗分为2行1列
stem（n,x,'.'）;%画脉冲图
title（'时域波形'）;
xlabel（'n'）;
ylabel（'x（n）'）;
subplot（2,1,2）;
stem（w/pi,mag）;%归一化角频率
axis（[0.20.502]）;%控制坐标范围以使谱线幅度合适
title（'1000点DFT频谱分析'）;
xlabel（'数字频率'）;
ylabel（'X（k）'）;
gridon;
实验结果：
2）DTMF信号频谱分析代码
clear
closeall%关闭所有窗口
fh=[1209,1336,1477,1633];%一个数组存放高频
fl=[697,770,852,941];%一个数组存放低频
fs=8000;%采样频率
N=1024;%采样点数
ts=1/fs;%采样周期
n=0:
N-1;
f=0:
fs/N:
fs/N*（N-1）;
key=zeros（10,N）;%生成0矩阵
key（1,:
）=cos（2*pi*fh
（1）*ts*n）+cos（2*pi*fl
（1）*ts*n）;%数字键1
key（2,:
）=cos（2*pi*fh
（2）*ts*n）+cos（2*pi*fl
（1）*ts*n）;%数字键2
key（3,:
）=cos（2*pi*fh（3）*ts*n）+cos（2*pi*fl
（1）*ts*n）;%数字键3
key（4,:
）=cos（2*pi*fh
（1）*ts*n）+cos（2*pi*fl
（2）*ts*n）;%数字键4
key（5,:
）=cos（2*pi*fh
（2）*ts*n）+cos（2*pi*fl
（2）*ts*n）;%数字键5
key（6,:
）=cos（2*pi*fh（3）*ts*n）+cos（2*pi*fl
（2）*ts*n）;%数字键6
key（7,:
）=cos（2*pi*fh
（1）*ts*n）+cos（2*pi*fl（3）*ts*n）;%数字键7
key（8,:
）=cos（2*pi*fh
（2）*ts*n）+cos（2*pi*fl（3）*ts*n）;%数字键8
key（9,:
）=cos（2*pi*fh（3）*ts*n）+cos（2*pi*fl（3）*ts*n）;%数字键9
key（10,:
）=cos（2*pi*fh
（2）*ts*n）+cos（2*pi*fl（4）*ts*n）;%数字键0
figure;
fori=1:
9%画出数字键1到9的频谱
subplot（3,4,i）
plot（f,abs（fft（key（i,:
））））;%利用fft求各个键的频谱
axis（[500170001000]）;%限定坐标轴范围
xlabel（'频率（hz）'）;ylabel（'幅度'）;%横纵坐标标识
title（i）;
grid;
subplot（3,4,10）%画出数字键0的频谱图
plot（f,abs（fft（key（10,:
））））;
axis（[500170001000]）;
xlabel（'频率（hz）'）;ylabel（‘幅度'）;
title（0）;
grid;end
DTMF信号频谱分析结果
实验二：
DTMF信号的编码
1、实验内容及要求
1）把您的联系电话号码通过DTMF编码生成为一个.wav文件。
³技术指标：
±根据ITUQ.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滤波器可以有效地提高计算效率。
其传递函数为：
³
2、实验目的
（a）复习和巩固IIR数字滤波器的基本概念；
（b）掌握IIR数字滤波器的设计方法；
（c）掌握IIR数字滤波器的实现结构；
（d）能够由滤波器的实现结构分析滤波器的性能（字长效应）；
（e）了解通信系统电话DTMF拨号的基本原理和IIR滤波器实现方法。
3.设计思路及实验步骤
1）编码：
DTMF 信号是将拨号盘上的 0~F 共16 个数字，用音频范围的 8 个频率来表示的一种编码方式。
8 个频率分为高频群和低频群两组，分别作为列频和行频。
每个字符的信号由来自列频和行频的两个频率的正弦信号叠加而成，组合方式如表格所示
1209hz
1336hz
1477hz
16333hz
697hz
1
2
3
A
770hz
4
5
6
B
852hz
7
8
9
C
941hz
*
0
#
D
根据表格即可得到各个数字对应的DTFM信号。
通过zeros全零矩阵来设置占空比，以达到题目要求。
得到信号后，使用sound函数来播放拨号音，writewave将信号写入声音文件。
2）解码：
在解码时，使用Goertzel算法。
滤波器调谐到这8个频率之上后，在相应的频率上的频谱值最大，通过与标准值的对比找出在DTMF图中的行和列，再对应出相应的拨号数字。
由Mock确定的Goertzel算法参数，假设一次谐波的DFT长度为205，此时对应的离散频率点k的值为18、20、22、24、31、34、38、42。
依次对应频率697hz,770hz,852hz,941hz,1209hz,1336hz,1447hz,1633hz。
4.实验代码及实验结果
1）实验代码：
（我的电话号码是）
N=400;fs=8000;%每个号码100ms
tm=[49,50,51,65;52,53,54,66;55,56,57,67;42,48,35,68];% DTMF表中键的16个ASCII码
n=1:
N;%取样点
fl=[697770852941];%低频
fh=[1209133614771633];%高频
x01=sin（2*pi*fl
（1）*n/fs）+sin（2*pi*fh
（1）*n/fs）;%1编码过程
x02=sin（2*pi*fl
（1）*n/fs）+sin（2*pi*fh（3）*n/fs）;%3
x03=sin（2*pi*fl
（2）*n/fs）+sin（2*pi*fh
（2）*n/fs）;%5
x04=sin（2*pi*fl
（1）*n/fs）+sin（2*pi*fh
（2）*n/fs）;%2
x05=sin（2*pi*fl
（1）*n/fs）+sin（2*pi*fh
（1）*n/fs）;%1
x06=sin（2*pi*fl（3）*n/fs）+sin（2*pi*fh（3）*n/fs）;%9
x07=sin（2*pi*fl
（1）*n/fs）+sin（2*pi*fh
（1）*n/fs）;%1
x08=sin（2*pi*fl
（2）*n/fs）+sin（2*pi*fh（3）*n/fs）;%6
x09=sin（2*pi*fl
（1）*n/fs）+sin（2*pi*fh（3）*n/fs）;%3
x10=sin（2*pi*fl（4）*n/fs）+sin（2*pi*fh
（2）*n/fs）;%0
x11=sin（2*pi*fl（3）*n/fs）+sin（2*pi*fh
（1）*n/fs）;%7
x=[x01,x02,x03,x04,x05,x06,x07,x08,x09,x10,x11];%组成矩阵
x01_z=[x01,zeros（1,400）];%补零满足占空比要求
x02_z=[x02,zeros（1,400）];
x03_z=[x03,zeros（1,400）];
x04_z=[x04,zeros（1,400）];
x05_z=[x05,zeros（1,400）];%补零满足占空比要求
x06_z=[x06,zeros（1,400）];
x07_z=[x07,zeros（1,400）];%补零满足占空比要求
x08_z=[x08,zeros（1,400）];
x09_z=[x09,zeros（1,400）];
x10_z=[x10,zeros（1,400）];
x11_z=[x11,zeros（1,400）];
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）;%发出声音
'phone.wav';%设置文件名
audiowrite（）;%写入声音文件
k=[1820222431343842];%解码
N=205;%取依次谐波的DFT点数
subplot（2,1,2）;
disp（['解码得到的号码是：
']）%在命令行窗口显示
fori=1:
11%对11位的已经编好码的电话号码进行解码
m=400*（i-1）;
X=goertzel（x（m+1:
m+N）,k+1）;%goertzel算法做变换
v=abs（X）;%求模
stem（k,v,'.'）;%画脉冲图
grid;
xlabel（'k'）;
ylabel（'x（k）'）;%作出8个点处X（K）的模值，比较最大的两个组合得到的数字即解码结果。
set（gcf,'color','w'）;%设置图片显示背景为白色
shg;%显示图形窗口
pause;
limit=80;
fors=5:
8;
ifv（s）>limit,break,end%查找列号
end
forr=1:
4;
ifv（r）>limit,break,end%查找行号
end
disp（[setstr（tm（r,s-4））]）%显示解码的字符
end
2）实验结果
K值为18和31，对应数字1
K值为18和38，对应数字3
K值为20和34，对应数字5
K值为18和34，对应数字2
K值为18和31，对应数字1
K值为22和38，对应数字9
K值为18和31，对应数字1
K值为20和38，对应数字6
K值为18和38，对应数字3
K值为24和34，对应数字0
K值为22和31，对应数字7
在命令窗口显示对应号码
实验三：
FIR数字滤波器的设计和实现
1、实验内容及要求：
³录制自己的一段声音，长度为10秒，取样频率32kHz，然后叠加一个高斯白噪声，使得信噪比为20dB。
请采用窗口法设计一个FIR带通滤波器，滤除噪声提高质量。
⏹提示：
³滤波器指标参考：
通带边缘频率为4kHz，阻带边缘频率为4.5kHz，阻带衰减大于50dB；
³Matlab函数y=awgn（x,snr,'measured'），首先测量输入信号x的功率，然后对其叠加高斯白噪声；
2、实验目的
³通过本次实验，掌握以下知识：
±FIR数字滤波器窗口设计法的原理和设计步骤；
±Gibbs效应发生的原因和影响；
±不同类型的窗函数对滤波效果的影响，以及窗函数和长度N的选择。
3.设计思路及实验步骤
1）fir滤波器设计原理：
设计fir滤波器的方法有很多，窗函数设计法，频率采样法，最优化设计法等，该实验我们采用窗函数设计法，其基本原理是用一定宽度的窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列。
其设计步骤主要有
·通过傅里叶逆变换获得理想滤波器的单位脉冲响应。
·由性能指标确定窗函数的种类和窗口长度
·求得实际滤波器的单位脉冲响应h（n）,即为所设计的FIR滤波器系数向量
·检验滤波器性能。
2）窗函数的选取：
滤波器的指标要求为 通带边缘为4KHZ，阻带边缘频率为4.5Khz，阻带衰减大于50db。
接下来应选取合适的窗函数实现滤波指标，然后计算出单位脉冲响应，再利用函数freqz（）求出其幅频特性看是否满足滤波器要求。
由于直接对信号加矩形窗截断会产生频谱泄露，为了改善频率泄露的问题，通常加非矩形窗，一般加汉明窗，因为汉明窗的幅频特性是旁瓣衰减较大，主瓣峰值与第一个旁瓣衰减可达50db。
在实验中使用Fir1（）设计满足滤波器指标要求的滤波器。
3）其次需要实现在声音上面叠加噪声，使用awgn函数。
使用作图函数分别画出原始声音信号和加噪声后的声音信号的时域和频域波形。
4）滤波操作：
使用matlab库函数hamming来进行滤波。
再作出滤波后的图形。
最后作出滤波器的幅频特性图。
4.实验代码及实验结果
1）实验代码
closeall;%关闭所有图形窗口
%原始声音信号频谱
[s1,fs]=audioread（'ashi.wav'）;%在当前文件夹下读取音频文件
fs=32000;%采样频率32khz
N=4096;%fft点数
k=1:
4096;%一个K对应一个数字频率
s1k=fft（s1,N）;%对原始音频信号进行频谱分析
figure
（1）;%在图形文件1中显示原始音频信号的时域和频谱
subplot（2,1,1）;
plot（s1）;
title（'原始信号（时域'）;
subplot（2,1,2）;
plot（32/4096*k,abs（s1k））;%频率分辨率=fs/N
axis（[0607]）;%设定横纵坐标范围
xlabel（'f（kHz）'）;
ylabel（'s1（k）'）;
title（'原声音信号（频域）'）;
%叠加噪声
s2=awgn（s1,20,'measured','db'）;%计算功率并叠加噪声
s2k=fft（s2,N）;%对叠加噪声过后的信号进行fft变换
figure
（2）;%在图形文件2中绘制叠加噪音后的时域频域幅频特性
subplot（2,1,1）;
plot（s2）;
title（'时域（叠加噪声）'）;
'noise.wav';
audiowrite（）;%将噪声文件写入
subplot（2,1,2）;
plot（32/4096*k,abs（s2k））;%频率分辨率=fs/N
title（'频域（叠加噪声后）'）;
xlabel（'f（kHz）'）;
ylabel（'s2（k）'）;
axis（[0607]）;
%设计滤波器过程
fp=4000;%通带起始频率
fr=4500;%阻带频率
wp=2*pi*fp/fs;%将滤波器指标转换到数字域
wr=2*pi*fr/fs;
width=wr-wp;%带宽
N1=ceil（6.6*pi/width）+1;%取整数
n=0:
1:
N1;
wc=（wr+wp）/2;%计算wc
alpha=（N1-1）/2;
n=0:
1:
N1-1;
m=n-alpha+eps;
hd=sin（wc*m）./（pi*m）;
b=fir1（N1,wc/（4*pi/3.65））;%返回截止频率为wc的N1阶滤波器hn，汉明窗
%通过Fir函数设计出满足要求的滤波器
w_hamming=（hamming（N1））';
h=hd.*w_hamming;%加窗，即时域进行卷积
[H,w]=freqz（h,[1],1000,'whole'）;%求幅频响应（fir滤波器的系统函数分子项为h，分母项为1）
H=（H（1:
1:
501））';
w=（w（1:
1:
501））';
mag=abs（H）;%求绝对值幅频响应
db=20*log10（（mag+eps）/max（mag））;%采用衰减特性，db做单位
pha=angle（H）;%相位响应
delta_w=2*pi/1000;
s3=filter（b,1,s2）;%对s2进行滤波，s3为滤波后的信号
'ashifil.wav';
audiowrite（）;%将滤波后的声音文件写入
%进行滤波处理
figure（3）;在图形文件3中绘制滤波后的波形文件的时域与频域响应
subplot（2,1,1）;
plot（s3）;
title（'时域（滤波后）'）;
s3k=fft（s3,N）;%对滤波后的信号进行频谱分析
subplot（2,1,2）;
plot（32/N*k,abs（s3k））;%分辨率
axis（[0607]）;
xlabel（'f（kHz）'）;
ylabel（'s3（k）'）;
title（'频域（滤波后）'）;
%滤波器特性
figure（4）;%在图形文件4中绘制滤波器的幅频特性
plot（w/pi,db）;
axis（[01-1500]）;
title（'滤波器（db）'）;
xlabel（'f'）;
ylabel（'deg'）;
2）实验结果截图
·原声音信号（时域和频域）
·叠加噪声后的声音信号（时域和频域）
·滤波后时域与频域信号
·滤波器幅频特性曲线
四．总结与感受
数字信号处理软件实验通过将课本上的知识与软件相结和，主要使用matlab作为工具，将理论课程进行了拓展和延伸。
第一个实验让我们通过matlab内置fft函数对数字信号进行频谱分析，在此过程中，最重要的是fft点数的选取及窗函数幅值的设定，该实验让我们掌握了用傅立叶变换进行信号分析时基本参数的选择。
了解了（DTFT）和（DFT）的区别，让我们掌握了离散傅立叶变换的基本原理、特性，体会了FFT的效率。
第二个实验通过将自己的电话号码进行编码和解码，让我们了解了DTMF信号的产生方式，以及熟悉了采用Goertzel算法解码，也掌握了采用该算法与fft算法解码的优劣，第三个实验通过将自己录制的语音信号加入白噪声然后观察其幅频响应并通过FIR滤波器滤除白噪声再观察其幅频响应，让我们掌握了FIR数字滤波器窗口设计法的原理和设计步骤,并体会到窗函数的类型和长度对结果的影响。
实验过程中遇到了许多问题，一方面是对matlab的不熟悉导致的编程难度大，另一方面是由于理论知识不牢固而导致的设计程序难度大。
这次实验中我对一些理论上的知识开始有了更深刻的认知，以前好像只是课本上的文字，读过也没有太大映像，而这次软件操作让我更深刻的理解了分辨率以及窗函数这两个知识点，也许还有更多纠正我以前错误的思维或是导通我迷糊的思维。
相信会对今后的专业学习和实践起到很大的帮助。
升级会员 