dsp-软件实验报告(精)

数字信号处理

MATLAB仿真实验报告

学院：电子工程学院

班级： 2011211203

学号： 2011210876

姓名：孙月鹏

班内序号： 04

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

、实验内容及要求

(1 离散信号的频谱分析：

设信号

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

(2 DTMF 信号频谱分析

用计算机声卡采用一段通信系统中电话双音多频（DTMF）拨号数字 0～9的数据，采用快速傅立叶变换（FFT）分析这10个号码DTMF拨号时的频谱。

2、实验结果

x(n的时域图与频谱：得到三根清晰的谱线

号码9的频谱号码8的频谱

号码7的频谱号码6的频谱

、实现代码及分析

(1第一小题：

k=1000; %DFT点数

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; %进行K点DFT变换

w=2*pi/k*[0:1:k-1]; %数字角频率

subplot(2,1,2;stem(w/pi,abs(xk; %画出频谱图

axis([0.2,0.5,0,2]; %设置窗函数的宽度与限幅

title('1000点dft';xlabel('数字频率';ylabel('|xk(k|';

% 此题关键在于DFT点数N的确定。经过计算和实验，当N=1000时能满足题目要求，看到3条清晰地谱线

（2）第二小题

clear;

close all;

f=[941 1336;697 1209;697 1336;697 1477;

770 1209;770 1336;770 1477;852 1209;

852 1336;852 1477] %0-9的频率

n=1:400;fs=4000; %取样频率为4000hz

fprintf('请输入数字(0 to 9：\n'

k=input (''

f1=f(k+1,1; %因为从0开始计算，+1得输

f2=f(k+1,2; %入数字的两个频率

N=400;

x1=sin(2*pi*f1*n/fs+sin(2*pi*f2*n/fs;%DTMF的输入信号时域

xn=[x1,zeros(1,400]; %补零

subplot(2,1,1;plot(xn %画出时域图

xlabel('n'

ylabel('xn'

subplot(2,1,2;

fn=fs*n/N; %取样点的频率

plot(fn,abs(fft(xn(1:400; %400点fft变换，画出频谱图

axis([0,4000,0,300]

xlabel('f'

ylabel('FFT'

二、实验二： DTMF 信号的编码

1、实验内容及要求

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

技术指标：

根据 ITU Q.23 建议，DTMF 信号的技术指标是：传送/接收率为每秒 10 个号码，或每个号码 100ms。

每个号码传送过程中，信号存在时间至少 45ms，且不多于 55ms，100ms 的其余时间是静音。

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

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

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

2、实验结果

输入号码界面

生成的时域图（占空比为50%）

电话号码的FFT图：每个数字都有两条主谱线

显示检测到的号码

3）实验代码及分析

%第一部分，产生编码

clc;

tm=[49,50,51,65; 52,53,54,66;55,56,57,67;42,48,35,68]; %DTMF表中的ASCII码

f1=[697,770,852,941]; %行频率向量

f2=[1209,1336,1477,1633]; %列频率向量

d=input('please enter number: ','s' %输入电话号码

sum=length(d; %电话号码长度

total_x=[]; %电话号码信号

sum_x=[];

sum_x=[sum_x,zeros(1,800];

for a=1:sum %循环sum次

symbol=abs(d(a; %求输入的ASCII码

for p=1:4;

for q=1:4;

if tm(p,q==abs(d(a; break,end %检测码相符的列号q

end

if tm(p,q==abs(d(a;break,end %检测码相符的行号p

end

n=1:400;

x=sin(2*pi*n*f1(p/8000+sin(2*pi*n*f2(q/8000;%构成双频信号x=[x,zeros(1,400]; %加长序列，增加静音

sum_x=sum_x+x;

total_x=[total_x,x]; %将所有编码连接起来

end

sound(total_x; %播放声音

t=total_x/2;

wavwrite(t,'我的手机号码'; %生成声音文件

plot(total_x;

title('DTMF信号时域波形';

%代码主要分成三部分，即根据输入的数字确定双频率、产生正弦信号和生成文件以及绘图。在产生正弦信号的过程中有增加静音和连接运算。整体上用循环对电话号码的每一位进行相同处理。

%第二部分，检测端

k=[18,20,22,25,32,35,38];

N=210;

tm=[49,50,51;52,53,54;55,56,57;0,48,0];

for i=1:sum j=800*(i-1;

X=goertzel(total_x(j+1:j+N,k+1; %算法公式

value=abs(X;

figure(2

subplot(2,6,i;

stem(k,value,'.','r'; %画出FFT频谱

title('FFT x(n';xlabel('k';ylabel('|X(k|';

limit=20;

for i1=5:7 %判决，确定频率

if value(i1>limit break;end

end

for j1=1:4

if value(j1>limit break;end

end

buffer(i=tm(j1,i1-4; %根据频率确定号码

end

disp(['接收端检测到的号码']

disp(setstr(buffer %显示检测到的号码

%检测端编码首先要根据goertzel公式算出FFT,再利用判决的方法找出频率的分布，从而确定输入的电话号码。

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

1、实验内容及要求：

录制自己的一段声音，长度为 45秒，取样频率 32kHz，然后叠加一个高斯白噪声，使得信噪比为 20dB。请采用窗口法设计一个 FIR 带通滤波器，滤除噪声提高质量。

? 提示：

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

Matlab 函数 y = awgn(x,snr,'measured' ，首先测量输入信号 x 的功率，然后对其叠加高斯白噪声；

2、实验结果

原始声音时域图与频谱图

叠加白噪声的时域图与频谱图

（但明显声音信号太弱，被白噪声盖过去了）

滤波器的幅频特性

滤波后的时域图与频谱特性

（可以看到45Khz以后的频谱都被滤掉了）

3、代码及分析

%产生原始声音及频谱

clear all

close all

[y1,fs,bits]=wavread('record-1.wav'; %按照采样频率取得y1 y2=y1(:,1;

fs=32000;

k=1:4096;

yk1=fft(y2,4096;

figure (1

subplot(2,1,1;plot(y2;title('原始声音时域';

subplot(2,1,2;plot(32/4096*k,abs(yk1;axis([0,17,0,0.05]; xlabel('f/kHZ';title('原始声音频域';

wavplay(y1,42000

%加噪及其频谱

y3=awgn(y2,20,'measured','db';%awgn?为加性高斯白噪声yk2=fft(y3,4096;

figure (2

subplot(2,1,1;plot(y3;title('加噪后声音时域';

subplot(2,1,2;plot(32/4096*k,abs(yk2;

xlabel('f/KHZ';title('加噪后声音频域';

wavplay(y3,42000;

%滤波器设计

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,y3;

figure(3

subplot(1,1,1;plot(w/pi,db;

title('滤波器幅频响应';

x=conv(h,y3;

%滤除噪声后频谱

figure(4

subplot(2,1,1;plot(x;

title('滤除噪声后时域图';

Xk=fft(x,4096;

subplot(2,1,2;

plot(32/4096*k,abs(Xk;axis([0,16,0,5]

title('滤除噪声后频域图'

xlabel('f/KHZ';

wavplay(x,42000;

%根据要求，所求的滤波器是一个低通滤波器

%因为衰减大于50DB，窗函数选择了汉明窗

4、错误分析与总结

问题主要出在第三道题目上：

1）实验要求的声音文件格式为wav，但录音软件的的格式wma，所以需要格式转换。首先实验的是直接更改格式名，尽管在播放器中可以播放，但仍然不可以识别；其次实验的是用matlab直接录音，但频谱集中在录音的采样频率附近，仍不能做正常的分析；最后用的是格式工厂的转换，尽管感觉频谱仍有改变，但并不影响实验结果。

2）用matlab对声音进行波形提取后，再用相同的fs播放，会发现声音的频率变低，进过实验，必须用比fs高一些的频率播放，才能听到正确的声音。

3）对比其他同学的频谱，发现我的原始声音的频谱的幅度就很低，加入白噪声之后，就很难看到声音的频谱了。但声音的播放就没有什么问题。直到最后也没有想明白为什么。

5、心得与体会

在这次实验中，我们分别用MATLAB实现了简单的FFT变换、IIR滤波器设计与FIR滤波器设计，真切的感受到了MATLAB仿真的便捷性。在MATLAB中，利用简单的几行代码，就能够完成特定的滤波器设计，是一种很好的辅助手段。于此同时，我也感受到了实践环节对知识的掌握、理论的学习非常重要。利用matlab

我可以很好的观察到N点的多少、窗函数的选择给FFT频谱、以及滤波器性能带来的改变。

随着课程的逐渐深入，我们大三学生将逐渐从课本上学习知识到与实践相结合。MATLAB是一项功能强大的软件，学完这门实验之后，我希望自己可以抽出实践继续学习MATLAB，将它运用到更多学科、更多学习上面。