[整理]dsp-软件实验报告

--------------------------数字信号处理MATLAB仿真实验报告学院：电子工程学院班级：2011211203学号：2011210876姓名：孙月鹏班内序号：04一、实验一：数字信号的FFT分析--------------------------1、实验内容及要求(1)离散信号的频谱分析：设信号此信号的0.3pi和0.302pi两根谱线相距很近，谱线0.45pi的幅度很小，请选择合适的序列长度N和窗函数，用DFT分析其频谱，要求得到清楚的三根谱线。(2)DTMF信号频谱分析用计算机声卡采用一段通信系统中电话双音多频（DTMF）拨号数字0～9的数据，采用快速傅立叶变换（FFT）分析这10个号码DTMF拨号时的频谱。2、实验结果x(n)的时域图与频谱：得到三根清晰的谱线号码9的频谱号码8的频谱00010450303024().*cos(.)sin(.)cos(.)xnnnn--------------------------号码7的频谱号码6的频谱3、实现代码及分析(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;closeall;f=[9411336;6971209;6971336;6971477;7701209;7701336;7701477;8521209;8521336;8521477]%0-9的频率n=1:400;fs=4000;%取样频率为4000hzfprintf('请输入数字(0to9)：\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文件。技术指标：根据ITUQ.23建议，DTMF信号的技术指标是：传送/接收率为每秒10个号码，或每个号码100ms。每个号码传送过程中，信号存在时间至少45ms，且不多于55ms，100ms的其余时间是静音。2）对所生成的DTMF文件进行解码。由于只需要知道8个特定点的频谱值，因此采用一种称为Goertzel算法的IIR滤波器可以有效地提高计算效率。其传输函数为2、实验结果输入号码界面2/1121()12cos(2/)jkNkezHzkNzz--------------------------生成的时域图（占空比为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('pleaseenternumber:','s')%输入电话号码sum=length(d);%电话号码长度total_x=[];%电话号码信号sum_x=[];sum_x=[sum_x,zeros(1,800)];fora=1:sum%循环sum次symbol=abs(d(a));%求输入的ASCII码forp=1:4;forq=1:4;iftm(p,q)==abs(d(a));break,end%检测码相符的列号qendiftm(p,q)==abs(d(a));break,end%检测码相符的行号pendn=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];%将所有编码连接起来endsound(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];fori=1:sumj=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('FFTx(n)');xlabel('k');ylabel('|X(k)|');limit=20;fori1=5:7%判决，确定频率ifvalue(i1)limitbreak;endendforj1=1:4ifvalue(j1)limitbreak;endendbuffer(i)=tm(j1,i1-4);%根据频率确定号码enddisp(['接收端检测到的号码'])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、代码及分析%产生原始声音及频谱clearallcloseall[y1,fs,bits]=wavread('record-1.wav');%按照采样频率取得y1y2=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，窗函数选择了汉明窗--------------------------四、错误分析与总结问题主要出在第三道题目上：1）实验要求的声音文件格式为wav，但录音软件的的格式wma，所以需要格式转换。首先实验的是直接更改格式名，尽管在播放器中可以播放，但仍然不可以识别；其次实验的是用matlab直接录音，但频谱集中在录音的采样频率附近，仍不能做正常的分析；最后用的是格式工厂的转换，尽管感觉频谱仍有改变，但并不影响实验结果。2）用matlab对声音进行波形提取后，再用相同的fs播放，会发现声音的频率变低，进过实验，必须用比fs高一些的频率播放，才能听到正确的声音。3）对比其他同学的频谱，发现我的原始声音的频谱的幅度就很低，加入白噪声之后，就很难看到声音的频谱了。但声音的播放就没有什么问题。直到最后也没有想明白为什么。五、心得与体会在这次实验中，我们分别用MATLAB实现了简单的FFT变换、IIR滤波器设计与FIR滤波器设计，真切的感受到了MATLAB仿真的便捷性。在MATLAB中，利用简单的几行代码，就能够完成特定的滤波器设计，是一种很好的辅助手