青岛科技大学《数字信号处理》课程设计(一)

1、青 岛 科 技 大 学数字信号处理课程设计报告数字信号分析及数字滤波器设计题 目 _张淑军指导教师_学生姓名_学生学号_院（部）_专业班_2014_年 _12_月 _25_日1. 目的与要求1.进一步巩固数字信号处理中的基本原理与方法，提高分析、解决实际问题的能力。 2.熟练掌握一门计算机语言，进行数字信号处理应用的开发设计，训练基本技能，如查阅设计资料和手册、程序的设计、调试等。3选择的题目为（一）数字信号分析及数字滤波器设计1.用以下方式产生三个不同频段的信号：（1）自己录制一段正常的语音文件；（2）录制一段环境噪声文件；（3）利用MATLAB产生一个不同于以上频段的信号。2.对上述三个信

2、号，进行频谱分析，画出三路信号的时域波形和频谱图，对进行对比分析。3.根据三路信号的频谱特点得到性能指标，由性能指标设计三个滤波器，并画出各滤波器的频域响应。4.将三路信号叠加为一路信号。5.用自己设计的滤波器对合成的信号进行滤波，分析得到信号的频谱，并画出滤波后信号的时域波形和频谱。2.主要技术和原理2.1语音的输入、处理以及信号的叠加 Matlab中语音信号的输入为x,fs,bits=wavread('C:UserslenovoDesktop原始语音.wav',50000 60000);其中wavread函数读取Microsoft的wav格式音频文件，x为自己随意定义的一个

3、输出函数，fs 是采样频率 bits 表示每个样点的位数，50000 60000表示读取从50000到60000个点到x中。X_voice=x(:,1)因为音频信号都有两个声道，所以提取出来一个声道进行分析Plot（x）为画出当前的波形，freqz（x）是求离散系统频响特性的函数。随机产生的信号（2）信号的叠加本次实验需要将正常语音信号与噪音信号还有matlab产生的一个信号进行叠加，生成一个新的信号，代码如下f1=22000; f2=25000; fs=50000; t=1/fs; length_x=length(x); input_x=zeros(1,length_x);for i=1:1

4、:length_x input_x(i)=sin(2*pi*f1*i*t)+sin(2*pi*f2*i*t);input_x=input_x'end 正常语音信号的代码：x,fs,bits=wavread('C:UserslenovoDesktop原始语音.wav',50000 60000);X_voice=x(:,1);input_x=X_voice;噪声语音的代码：x1,fs1,bits1=wavread('C:UserslenovoDesktopÔëÒô1.wav',50000 60000);X3_voice

5、=x1(:,1);input_x1=X3_voice;但是信号叠加有一个非常重要的条件就是点数必须相同，不仅点数相同而且点数矩阵得是同一类型，在原始语音和噪声语音的点数10001个，但他们产生的是10001行1列矩阵，而matlab产生的正弦信号却是1行10001列矩阵，结果也不能叠加，这就需要将1行10001列矩阵变成10001行1列矩阵，input_x=input_x'反转语句就是讲input_x的行和列互换。叠加语句：for i=1:1:length_xy1=X_voice+X3_voice+input_x;end2.2FIR滤波器设计原理时间窗口设计法是从单位脉冲响应序列着手，

6、使h(n)逼近理想的单位脉冲响应序列hd(n)。我们知道hd(n)可以从理想频响通过付氏反变换获得但一般来说，理想频响是分段恒定，在边界频率处有突变点，所以，这样得到的理想单位脉冲响应hd(n)往往都是无限长序列，而且是非因果的。但FIR的h(n)是有限长的，问题是怎样用一个有限长的序列去近似无限长的。最简单的办法是直接截取一段hd(n)代替 h(n) 。这种截取可以形象地想象为h(n)是通过一个“窗口”所看到的一段hd(n)，因此 ，h(n)也可表达为h(n)和一个“窗函数”的乘积，即 h(n)=w(n) hd(n) 在这里窗口函数就是矩形脉冲函数RN（n），当然以后我们还可看到，为了改善设

7、计滤波器的特性，窗函数还可以有其它的形式，相当于在矩形窗内对hd(n)作一定的加权处理。（2）本次实验只采用了矩形窗设计方法，由于hd(n)是无限长的，故要用一个有限长的“窗函数”序列w(n)将hd(n)加以截断，窗的点数是N点。截断后的序列为hd(n) h(n)=w(n) hd(n)， 0<=n<=N-1窗的点数N及窗的形状是两个极其重要的参数。加窗后实际的频率响应 =DTFTh(n)=DTFT=(3)FIR带通、高通、带阻、多通带多阻带滤波器的设计（1）N长的FIR全通滤波器h（n）的函数表达当，则低通滤波器变为全通滤波器： （2）N长的FIR高通滤波器h（n）的函数表达式由高

8、通的频谱结构可知，高通滤波器的频谱=全通滤波器的频谱-低通滤波器的频谱。（3）N长的FIR带通滤波器h（n）的函数表达式带通滤波器的频谱=低通滤波器的频谱-低通滤波器的频谱。（3）N长的FIR带通滤波器h（n）的函数表达式带阻滤波器的频谱=全通滤波器的频谱-带通滤波器的频谱。 （4）N长的FIR带阻滤波器h（n）的函数表达式 带阻滤波器的频谱=全通滤波器的频谱-带通滤波器的频谱 （5）N长的FIR多通带多带阻滤波器h（n）的函数表达式 多通带多阻带滤波器的频谱=全通滤波器的频谱-带通滤波器的频谱-带通滤波器的频谱。3.总体方案设计3.1开发环境和工具 Pc window8.1 matlab20

9、12a3.2具体框架Matlab产生一个信号噪声语音的录入原始信号的录入绘制时域频域波形绘制时域频域波形绘制时域频域波形三路信号叠加成一路信号绘制时域频域波形设计多通多阻滤波器设计高通滤波器设计低通滤波器与y1相卷积进行滤波与y1相卷积进行滤波与y1相卷积进行滤波语音恢复成wav格式语音恢复成wav格式语音恢复成wav格式3.3设计过程 （1）信号的获取采集（原始语音，噪声，matlab产生的信号） 录取一段语音信号，fs=44100赫兹，共抽样10001个点录取一段噪声语音，本次噪声语音就是用电脑自己录的杂音当做噪声，且fs=44100，也是抽取了10001个点用matlab产生了一个正弦信

10、号抽样频率为50000，f1和f2分别为22000和25000赫兹。 （2）对做好的信号要进行处理，用wavread函数将语音信号读进来，取他们的一列，对于正弦信号要进行转置input_x=input_x'以便于后面对信号的叠加做准备。 （3）绘制语音信号的时域波形的频域波形 （4）根据（3）设计出来的频谱特性进行分析，根据分析的结果设计滤波器，在画出各个滤波器的频谱以及时域波形。 （5）将做好的三路信号叠加成一路信号，在对叠加好的信号画频谱。 （6）将设计好的滤波器分别与此信号卷积进行滤波，在画出滤波后的信号的频谱。 （7）最后将语音恢复，输出到指定文件夹下。4.实验结果与分析4.1

11、这部分是实现对原始语音的处理具体代码如下：x,fs,bits=wavread('C:UserslenovoDesktop原始话音.wav',50000 60000); %对原始话音读取到x中X_voice=x(:,1);%提取声道1input_x=X_voice;%将声道1赋给input_xfigure(1);plot(input_x);对input_x画时域波形title('原始话音波形');figure(1);freqz(input_x);%对input_x画频谱波形title('原始话音频谱');X1_voice=x(:,2);%提取声道2

12、input_x=X1_voice;%将声道2赋给input_xfigure(3);plot(input_x);title('原始话音波形');figure(4);freqz(input_x);title('原始话音频谱');y=X_voice+X1_voice;&将两个信号叠加为原始语音信号figure(5);plot(y);title('原始话音波形');figure(6);freqz(y);title('原始话音频谱');主要用到的函数wavread是对语音信号的读取，但只能是wav格式，plot(input_x);是

13、画图函数，将input_x画出来，freqz(input_x);是求input_x的频谱并将其画出。实验结果图如下：上图为声道1的波形和频谱，只是提取了10001个点进行分析，此频谱特点是图上在f=2800,11250,21400赫兹处出现了尖峰。这里只对提取的一个声道分析，代表了此信号。4.2实现对噪音信号的处理模块及其实现代码：x1,fs1,bits1=wavread('C:UserslenovoDesktopÔëÒô1.wav',50000 60000);X3_voice=x1(:,1);input_x1=X3_voice;figu

14、re(9);plot(input_x1);title('原始噪音波形');figure(10);freqz(input_x1);title('原始噪音频谱图');此模块与正常语音的的处理大同小异，只是在噪音里我只提取了一个声道，用到的函数与处理正常语音用到的函数一样。实验结果图如下以下两个图一个是噪音频谱，另一个是噪音波形，我也是对噪音提取了10001个点进行分析，从频谱上看，在f=6700赫兹以前信号的频率比较高，但不是很高，就设计一个高通滤波器将较高的信号滤除，使产生一个平稳的信号，此噪音不突出，一方面跟采样点数有关系，另一方面录制噪音的时候尖峰不够，但主要

15、考察滤波器设计及其分析，不影响整体。4.3实现matlab产生的正弦信号及其代码：f1=22000; % 输入信号f1f2=25000; % 输入信号f2fs=50000; % fs为抽样频率 t=1/fs; length_x=length(x); %输入信号的长度input_x=zeros(1,length_x);for i=1:1:length_x input_x(i)=sin(2*pi*f1*i*t)+sin(2*pi*f2*i*t);input_x=input_x'%将input_x转置end figure(7); plot(input_x); title('产生信号的

16、时域波形');grid on;figure(8); freqz(input_x); title('产生信号的频谱波形');此模块先设置了两个信号的频率，抽样频率要大于等于2倍的f1和f2，i=1:1:length_x是从1开始到length_x逐步取值赋给x，且步长为1，input_x=input_x'是对input_x进行转置本来是100011矩阵转置后变为110001矩阵，grid on是画图时出现网格，其他函数以上均有说明。实验结果图如下：以上两个图形是产生信号的波形及频谱，由于采样点数很大，导致时域波形看起来就是一个蓝色矩形，点数少了能看见是一个正弦信号

17、。从频谱上，在f=21000赫兹处出现了一个尖峰。4.4实现信号的叠加的代码：for i=1:1:length_x y1=X_voice+X3_voice+input_x;endfigure(11);freqz(y1);title('叠加后信号的频谱');信号叠加要考虑到长度相同，且矩阵也要相同。实验结果图如下：综合原始语音，噪声语音，还有产生的信号的频谱分析，在0.9出出现了一个尖峰，原来的尖峰依旧存在，叠加信号的频谱跟原始语音信号非常相似，产生的信号叠加后比原始信号频谱多了一个尖峰，但不明显，相反噪音信号无尖峰，对正常语音信号影响并不大，可能是点数没有取好，但只要有影响就行

18、。后期处理好就行。而时域波形可以看见有很多蜂刺，那些蜂刺可能是噪音的影响。4.5由原始语音设计的滤波器对叠加后的信号处理及实现其代码：f_L1=2600;%其频率值f_H1=3000;f_L2=11000; f_H2=11500;f_L3=21900;f_H3=22600;w_L1=2*pi*f_L1/fs;%转化为数字频率值w_H1=2*pi*f_H1/fs;w_L2=2*pi*f_L2/fs; w_H2=2*pi*f_H2/fs;w_L3=2*pi*f_L3/fs;w_H3=2*pi*f_H3/fs;N=22001;M=N-1;h=zeros(1,N);for i=0:N-1%多通多阻滤波

19、器的设计 if i=M/2 h(i+1)=(pi-w_H1+w_L1-w_H2+w_L2-w_H3+w_L3)/pi; else h(i+1)=(sin(i-M/2)*pi)-sin(i-M/2)*w_H1)+sin(i-M/2)*w_L1)-sin(i-M/2)*w_H2)+sin(i-M/2)*w_L2)-sin(i-M/2)*w_H3)+sin(i-M/2)*w_L3)/(i-M/2)*pi); end end figure(16); plot(h);title('由原始语音设计的滤波器的频谱 ');figure(17);freqz(h); title('经滤波后

20、的y1');y3=conv(y1,h);freqz(y3);其中for循环是设计多通多阻，两个低通相减就是带通，再用全通减去带通就得到带阻，带通的代码sin(i-M/2)*w_H1)-sin(i-M/2)*w_L1) /(i-M/2)*pi)，全通的代码sin(i-M/2)*pi)/(i-M/2)*pi)，以上的频率是分析其尖峰所选择的频率。其滤波器频域波形及滤波后的时域波形、频域波形：上图是经过滤波后y1的频谱，从图上可以看出有尖峰的地方已经滤掉可以看出滤波器设计还算成功上图是设计的滤波器的频谱，含有多通多阻，从图中可以看出有三处是阻带，就是为了滤掉y1的三个尖峰，其他可以安全通过。

21、下图为滤波后y1的时域波形，波形大概成对称性4.5由噪音语音设计的滤波器对叠加后的信号处理及实现其代码：N=1001;M=N-1;wc=2*pi*6700/fs;%截至的角频率for i=0:N-1 %设计的高通滤波器 if i=M/2 h(i+1)=(pi-wc)/pi; % 当i=M/2时的极限值 else h(i+1)=(sin(i-M/2)*pi)-sin(i-M/2)*wc)/(i-M/2)*pi); % 高通滤波器的h(n)表达式 end end figure(12);freqz(h);title('由噪声设计的滤波器的频谱');y=conv(y1,h);figur

22、e(13);freqz(y);title('经过滤波后的y1')这是一个高通滤波器，高通滤波器=全通滤波器-低通滤波器，wc是滤波器的截止频率在此处为过渡段，其滤波器频域波形及滤波后的时域波形、频域波形：上图是设计的高通滤波器的频谱图上图是经过高通滤波器后y1的时域波形，高通对叠加信号的时域波形影响并不大只是尖峰变密了。上图是经过滤波后y1的频谱波形，可以看见将噪音信号的较高频率的滤掉了。4.5由matlab产生的信号设计的滤波器对叠加后的信号处理及实现其代码：N=1000;M=N-1;wc=2*pi*21000/fs;for i=0:N-1%设计的低通滤波器 if i=M/2

23、 h(i+1)=wc/pi; % 当i=M/2时的极限值 else h(i+1)=sin(i-M/2)*wc)/(i-M/2)*pi); % 低通滤波器的h(n)表达式end end figure(14);freqz(h); %显示低通滤波器的频率响应函数title('由产生的信号设计出的滤波器的频谱');y2=conv(y1,h);figure(15);freqz(y2);title('经过滤波后的y1');低通滤波器的函数表达为sin(i-M/2)*wc)/(i-M/2)*pi);其滤波器频域波形及滤波后的时域波形、频域波形：上图是根据产生信号频谱设计的滤波器的频谱此为低通滤波器。上图为滤波后的y1时域波形，此时域图跟叠加后时域波形变化巨大看上去更像正弦信号。上图是经过低通滤波器y1的频谱。可以看到将产生信号频谱的尖峰滤掉了。4.6