北邮数字信号处理MATLAB实验报告

1、数字信号处理软件实验MatLab仿真实验报告学院:电子工程学院班级：2013211202姓名：学号：实验一：数字信号的 FFT 分析1、实验内容及要求(1) 离散信号的频谱分析： 设信号 此信号的0.3pi 和 0.302pi两根谱线相距很近，谱线 0.45pi 的幅度很小，请选择合适的序列长度 N 和窗函数，用 DFT 分析其频谱，要求得到清楚的三根谱线。(2) DTMF 信号频谱分析 用计算机声卡采用一段通信系统中电话双音多频（DTMF）拨号数字 09的数据，采用快速傅立叶变换（FFT）分析这10个号码DTMF拨号时的频谱。 2、实验目的 通过本次实验，应该掌握：(a) 用傅立叶变换进行信

2、号分析时基本参数的选择。 (b) 经过离散时间傅立叶变换（DTFT）和有限长度离散傅立叶变换（DFT） 后信号频谱上的区别，前者 DTFT 时间域是离散信号，频率域还是连续的，而 DFT 在两个域中都是离散的。(c) 离散傅立叶变换的基本原理、特性，以及经典的快速算法（基2时间抽选法），体会快速算法的效率。(d) 获得一个高密度频谱和高分辨率频谱的概念和方法，建立频率分辨率和时间分辨率的概念，为将来进一步进行时频分析（例如小波）的学习和研究打下基础。(e) 建立 DFT 从整体上可看成是由窄带相邻滤波器组成的滤波器组的概念，此概念的一个典型应用是数字音频压缩中的分析滤波器，例如 DVD AC3

3、 和MPEG Audio。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拨号盘中的每一个按键均对应两个频率，一个高频

4、，一个低频，每个数字由两个不同频率的正弦波组成，高频为：1209HZ,1336HZ,1633HZ,1477HZ,低频为697HZ,770hz，852hz,941hz。分别用两个数组来装载高频频率和低频频率，然后再分别产生两路频率混叠的信号。再利用fft频谱分析做出频谱图。4.实验代码及实验结果1)离散信号的频谱分析代码clf;close all;%关闭所有图形窗口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/

5、N*0:1:N-1; %数字角频率wsubplot(2,1,1); %将图形窗分为2行1列stem(n,x,'.'); %画脉冲图title('时域波形');xlabel('n');ylabel('x(n)');subplot(2,1,2);stem(w/pi,mag); %归一化角频率axis(0.2 0.5 0 2); %控制坐标范围以使谱线幅度合适title('1000点DFT频谱分析');xlabel('数字频率');ylabel('X(k)');grid on;实验结果：2

6、) DTMF 信号频谱分析代码clearclose all%关闭所有窗口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);%数字键1key(2,:)=cos(2*pi*fh(2)*ts*n)+cos(2*pi*fl(1)*ts*n); %数字键2key

7、(3,:)=cos(2*pi*fh(3)*ts*n)+cos(2*pi*fl(1)*ts*n); %数字键3key(4,:)=cos(2*pi*fh(1)*ts*n)+cos(2*pi*fl(2)*ts*n); %数字键4key(5,:)=cos(2*pi*fh(2)*ts*n)+cos(2*pi*fl(2)*ts*n); %数字键5key(6,:)=cos(2*pi*fh(3)*ts*n)+cos(2*pi*fl(2)*ts*n); %数字键6key(7,:)=cos(2*pi*fh(1)*ts*n)+cos(2*pi*fl(3)*ts*n); %数字键7key(8,:)=cos(2*pi*

8、fh(2)*ts*n)+cos(2*pi*fl(3)*ts*n); %数字键8key(9,:)=cos(2*pi*fh(3)*ts*n)+cos(2*pi*fl(3)*ts*n); %数字键9key(10,:)=cos(2*pi*fh(2)*ts*n)+cos(2*pi*fl(4)*ts*n); %数字键0figure;for i=1:9%画出数字键1到9的频谱subplot(3,4,i)plot(f,abs(fft(key(i,:);%利用fft求各个键的频谱axis(500 1700 0 1000);%限定坐标轴范围xlabel('频率（hz）');ylabel('

9、幅度');%横纵坐标标识title(i);grid;subplot(3,4,10)%画出数字键0的频谱图plot(f,abs(fft(key(10,:);axis(500 1700 0 1000);xlabel('频率（hz）');ylabel(幅度');title(0);grid;endDTMF 信号频谱分析结果实验二： DTMF 信号的编码1、实验内容及要求1）把您的联系电话号码 通过DTMF 编码生成为一个 .wav 文件。³ 技术指标：± 根据 ITU Q.23 建议，DTMF 信号的技术指标是：传送/接收率为每秒 10 个号码，或每

10、个号码 100ms。± 每个号码传送过程中，信号存在时间至少 45ms，且不多于 55ms，100ms 的其余时间是静音。± 在每个频率点上允许有不超过 ±1.5% 的频率误差。任何超过给定频率 ±3.5% 的信号，均被认为是无效的，拒绝接收。（其中关键是不同频率的正弦波的产生。可以使用查表方式模拟产生两个不同频率的正弦波。正弦表的制定要保证合成信号的频率误差在±1.5%以内，同时使取样点数尽量少） 2）对所生成的DTMF文件进行解码。³ DTMF 信号解码可以采用 FFT 计算 N 点频率处的频谱值，然后估计出所拨号码。但 FFT计

11、算了许多不需要的值，计算量太大，而且为保证频率分辨率，FFT的点数较大，不利于实时实现。因此，FFT 不适合于 DTMF 信号解码的应用。³ 由于只需要知道 8 个特定点的频谱值，因此采用一种称为 Goertzel 算法的 IIR 滤波器可以有效地提高计算效率。其传递函数为：³2、实验目的(a)复习和巩固 IIR 数字滤波器的基本概念；(b)掌握 IIR 数字滤波器的设计方法；(c)掌握 IIR 数字滤波器的实现结构；(d)能够由滤波器的实现结构分析滤波器的性能（字长效应）；(e)了解通信系统电话 DTMF 拨号的基本原理和 IIR 滤波器实现方法。3.设计思路及实验步骤1

12、）编码：DTMF 信号是将拨号盘上的 0F 共16 个数字，用音频范围的 8 个频率来表示的一种编码方式。8 个频率分为高频群和低频群两组，分别作为列频和行频。 每个字符的信号由来自列频和行频的两个频率的正弦信号叠加而成，组合方式如表格所示1209hz1336hz1477hz16333hz697hz123A770hz456B852hz789C941hz*0#D根据表格即可得到各个数字对应的DTFM信号。通过zeros全零矩阵来设置占空比，以达到题目要求。得到信号后，使用sound函数来播放拨号音，writewave

13、将信号写入声音文件。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; %每个号码100mstm=49,

14、50,51,65;52,53,54,66;55,56,57,67;42,48,35,68; % DTMF表中键的16个ASCII码n=1:N;%取样点fl=697 770 852 941;%低频fh=1209 1336 1477 1633;%高频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); %3x03=sin(2*pi*fl(2)*n/fs)+sin(2*pi*fh(2)*n/fs); %5x04=sin(2*pi*fl(1)*n/

15、fs)+sin(2*pi*fh(2)*n/fs); %2x05=sin(2*pi*fl(1)*n/fs)+sin(2*pi*fh(1)*n/fs); %1x06=sin(2*pi*fl(3)*n/fs)+sin(2*pi*fh(3)*n/fs); %9x07=sin(2*pi*fl(1)*n/fs)+sin(2*pi*fh(1)*n/fs);%1x08=sin(2*pi*fl(2)*n/fs)+sin(2*pi*fh(3)*n/fs); %6x09=sin(2*pi*fl(1)*n/fs)+sin(2*pi*fh(3)*n/fs); %3x10=sin(2*pi*fl(4)*n/fs)+sin

16、(2*pi*fh(2)*n/fs); %0x11=sin(2*pi*fl(3)*n/fs)+sin(2*pi*fh(1)*n/fs); %7x=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,40

17、0); %补零满足占空比要求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);%发出声音filename='phone.wav'%设置文件名audiow

18、rite(filename,x_z,fs);%写入声音文件k=18 20 22 24 31 34 38 42;%解码N=205;%取依次谐波的DFT点数subplot(2,1,2);disp('解码得到的号码是：')%在命令行窗口显示for i=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)&

19、#39;);%作出8个点处X(K)的模值，比较最大的两个组合得到的数字即解码结果。 set(gcf,'color','w');%设置图片显示背景为白色 shg;%显示图形窗口 pause; limit=80; for s=5:8; if v(s)>limit,break,end%查找列号endfor r=1:4; if v(r)>limit,break,end%查找行号enddisp(setstr(tm(r,s-4)%显示解码的字符end2）实验结果 K值为18和31，对应数字1 K值为18和38，对应数字3 K值为20和34，对应数字5 K值为18

20、和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 带通滤波器，滤除噪声提高质量。n 提示：³ 滤波器指标参考：通带边缘频率为 4kHz，阻带边缘频率为4.5kHz，阻带衰减大于 50

21、dB；³ Matlab 函数 y = awgn(x,snr,'measured') ，首先测量输入信号 x 的功率，然后对其叠加高斯白噪声；2、实验目的³ 通过本次实验，掌握以下知识：± FIR 数字滤波器窗口设计法的原理和设计步骤；± Gibbs 效应发生的原因和影响；± 不同类型的窗函数对滤波效果的影响，以及窗函数和长度 N 的选择。3.设计思路及实验步骤1）fir滤波器设计原理：设计fir滤波器的方法有很多，窗函数设计法，频率采样法，最优化设计法等，该实验我们采用窗函数设计法，其基本原理是用一定宽度的窗函数截取无限脉冲响应

22、序列获得有限长的脉冲响应序列。其设计步骤主要有·通过傅里叶逆变换获得理想滤波器的单位脉冲响应。·由性能指标确定窗函数的种类和窗口长度·求得实际滤波器的单位脉冲响应h(n),即为所设计的FIR滤波器系数向量·检验滤波器性能。2）窗函数的选取：滤波器的指标要求为 通带边缘为4KHZ，阻带边缘频率为4.5Khz，阻带衰减大于50db。接下来应选取合适的窗函数实现滤波指标，然后计算出单位脉冲响应，再利用函数freqz()求出其幅频特性看是否满足滤波器要求。由于直接对信号加矩形窗截断会产生频谱泄露，为了改善频率泄露的问题，通常加非矩形窗，一般加汉明窗，因

23、为汉明窗的幅频特性是旁瓣衰减较大，主瓣峰值与第一个旁瓣衰减可达50db。在实验中使用Fir1()设计满足滤波器指标要求的滤波器。3）其次需要实现在声音上面叠加噪声，使用awgn函数。使用作图函数分别画出原始声音信号和加噪声后的声音信号的时域和频域波形。4）滤波操作：使用matlab库函数hamming来进行滤波。再作出滤波后的图形。最后作出滤波器的幅频特性图。4.实验代码及实验结果1）实验代码close all;%关闭所有图形窗口 %原始声音信号频谱 s1,fs=audioread('ashi.wav');%在当前文件夹下读取音频文件 fs=32000; %采样频率32khz

24、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(0 6 0 7);%设定横纵坐标范围 xlabel('f(kHz)'); ylabel('s1(k)'); title('原声音信号（频域）

25、'); %叠加噪声 s2=awgn(s1,20,'measured','db');%计算功率并叠加噪声 s2k=fft(s2,N);%对叠加噪声过后的信号进行fft变换 figure(2);%在图形文件2中绘制叠加噪音后的时域频域幅频特性 subplot(2,1,1); plot(s2); title('时域（叠加噪声）'); filename='noise.wav'audiowrite(filename,s2,48000);%将噪声文件写入 subplot(2,1,2); plot(32/4096*k,abs(s2k);

26、 %频率分辨率=fs/Ntitle('频域（叠加噪声后）');xlabel('f(kHz)');ylabel('s2(k)');axis(0 6 0 7);%设计滤波器过程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;%计算wcalpha=(N1-1)/2;n=0:1:N1-1;m=n-alpha+eps;hd=sin(wc*m

27、)./(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为滤波后的信号filename='ashifil.wav'audiowrite(filename,s3,48000);%将滤波后的声音文件写入%进行滤波处理figure(3);在图形文件3中绘制滤波后的波形文件的时域与频域响应subplot(2,1,1);plot(s3);title('时域（滤波后）');s3k=fft(s3,N);%对滤波后的信号进行频