青 岛 科 技 大 学

《数字信号处理》课程设计报告

题 目 ______________

 

 

 

 

指导教师 _{__________________________}

学生姓名 _{__________________________}

学生学号 _{__________________________}

 

_{_______________________________} 院（部） _{____________________________} 专业 _{________________} 班

_2014_ 年 _12_月 _21_日

1.目的与要求

数字信号处理课程设计

      现有原始语音文件“语音001.wav”，以及混有正弦噪声的语音文件“original_voice_Noise_by_SineWaves.wav”，试解决如下问题：

1）试分析有几个正弦噪声，他们的频率是多少Hz？

2）设计滤波器，一并滤除这些噪声，并将滤除噪声后的语音序列，保存为“voice_noise_filteredOUT.wav”文件，用播放器试听该文件，并与原始语音文件“语音001.wav”比较，感受噪声抑制的效果。

3）利用实验一的相关计算方法，大致判定这些正弦噪声在原始语音文件“语音001.wav”中出现的位置。

2.主要技术和原理

2.1  学会使用matlab软件，掌握利用matlab进行数字信号处理的方法。

2.2  学会如何利用matlab库函数进行声音文件的输入与输出存储。

2.3  如何利用matlab绘制各种图形。

2.4  无限长单位冲击响应（IIR）数字滤波器与有限长单位冲击响应（FIR）数字滤波器的设计。

2.5  快速傅里叶变换。

2.6  线性卷积。

2.7  幅频响应与相频响应的分析。

2.8  对滤波器滤波效果的检测。

2.9  数字信号处理方法是以数值计算为特点，它与模拟电路处理信号相比具有精度高、抗干扰能力强、便于集成、灵活性高等优点。

 

3.总体方案设计

3.1  声音文件的读入。

利用如下语句进行声音信号的读入：

    [x,fs,nbits] = wavread(‘D:\...\...\sound_and_noise3.wav’);

注意：声音文件路径要写对。

函数说明：（1）x为采样数据存放数组。在matlab的Command Window中通过输入whos x语句可以发现x为列向量，通过x = x’;将其转化为行向量。

         （2）fs为抽样频率，fs = 44100Hz。

         （3）nbits为采样数据位宽，nbits = 16bit。

3.2  时域观察。

利用plot（x）;语句可以画出采样信号的时域波形，但是由于采样点数太多,有207270个，所以在时域也观察不出有价值的信息。

3.3  频域分析。

利用freqz（x）函数，画出输入语音信号的幅频响应与相频响应，如下图：

观察幅频响应图，发现图中许多高峰，并不能定位噪音的频率，必须寻求其他的方法来确定噪音的频率，与同学交流后发现采用3.4中方法能准确的找出噪音的频率。

3.4    噪音频率的确定。

利用快速傅里叶变换，将时域信号x变换到频域分析。

     f = fs*(0:x_len/2-1)/x_len;        %模拟域频率 

X = fft(xdeal,x_len);              %快速傅里叶变换

plot(f,abs(X(1:x_len/2)));         %画图

   图中，三个最高峰即为噪音频率，原始信号中有3个不同频率的噪音，三个噪音的频率分别为：f1 = 3500Hz，f2 = 4600Hz，f3 = 5500Hz.

3.5 滤波器的设计

常用的滤波器有IIR（Infinite Impulse Response）数字滤波器与FIR（Finite Impulse Response filter）数字滤波器两种。

IIR滤波器在设计上可以借助成熟的模拟滤波器成果，如巴特沃斯（Butterworth）、契比雪夫（Chebyshev）、椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，现根据数字滤波器的指标，求出模拟滤波器的指标，然后设计出模拟滤波器，最后根据冲击响应不变法或双线性变换法，将模拟滤波器转化为数字滤波器。但是IIR数字滤波器的相位不好控制，对相位要求较高时，需加一个全通系统，对相位加以校正。

FIR滤波器采用窗函数法设计。与IIR滤波器相比，FIR滤波器有以下几个特点：

(1)由于单位抽样响应h(n)是有限长的，因此FIR数字滤波器一定是稳定的。

(2)经延时移位，h(n)总可变成因果序列，所以FIR DF（Digital Filter数字滤波器）总可以由因果系统实现。

(3)h(n)为有限长，可以用FFT实现FIR DF。

(4)FIR的系统函数是Z-1的多项式，是全零点系统，故IIR的设计方法不适用。

(5)FIR的相位特性可以是相性的，因此，它有更广泛的应用。

(6)FIR滤波器可以设计出理想正交变换器、理想微分器、线性调频器、多通带多阻带等各种网络。

综上所述，选用窗函数法设计FIR滤波器来去除噪声。矩形窗函数法设计各种滤波器的原理：

（1）  FIR低通滤波器：

（2）  FIR全通滤波器：

当，则低通变为全通： 

（3）  FIR高通滤波器：

        高通滤波器的频谱=全通滤波器的频谱-低通滤波器的频谱

 

（4）  FIR带通滤波器：

   带通滤波器的频谱=低通滤波器的频谱-低通滤波器的频谱。

（5）  FIR带阻滤波器：

       带阻滤波器的频谱=全通滤波器的频谱-带通滤波器的频谱

（6）  FIR多通带多阻带滤波器：

滤波器的频谱=全通滤波器的频谱-带通滤波器的频谱-带通滤波器的频谱

       

这里，需要设计一个四通带三阻带的FIR滤波器，阻带带宽设为100Hz，滤波器的阶数N=100，用matlab编程实现如下：

%%%%构建FIR滤波器

fL1 = 3450;        %噪声1的频率f1=3500Hz

fH1 = 3550;

fL2 = 4550;        %噪声2的频率f2=4600Hz

fH2 = 4650;

fL3 = 5450;        %噪音3的频率f3=5500Hz

fH3 = 5550;

wL1 = 2*pi*fL1/fs;    %化为数字角频率

wH1 = 2*pi*fH1/fs;

wL2 = 2*pi*fL2/fs;

wH2 = 2*pi*fH2/fs;

wL3 = 2*pi*fL3/fs;

wH3 = 2*pi*fH3/fs;

N = 1001;

M = (N-1)/2;

h = zeros(1,N);

for n=0:1:N-1

    if n==M;

        h(n+1) = (pi - wH1 + wL1 - wH2 + wL2 - wH3 + wL3)/pi;

    else

        h(n+1) = (sin(pi*(n-M)) - sin(wH1*(n-M)) + sin(wL1*(n-M)) - sin(wH2*(n-M)) + sin(wL2*(n-M)) - sin(wH3*(n-M)) + sin(wL3*(n-M)))/(pi*(n-M));

    end

end

利用plot(h)语句观察所构建滤波器的时域波形如下：

3.6 滤波器频域分析

滤波器的幅频响应表示系统的功能，即表示不同频谱成分的信号在该系统中的衰减程度，相频响应表示信号通过系统时，不同频谱成分的延时情况。相性相位条件下群延时的特点是群延时为常数，其作用是使系统通带内无失真传输的必要条件。

利用freqz（h）语句观察滤波器h(n)的幅频特性与相频特:

通过图形可以看出，该滤波器有三个阻带，且为线性相位（相频响应为一条直线）。

3.7 对信号滤波去噪

利用conv（）函数对x与h进行时域卷积，所得结果y_temp即为滤波后的信号。但这样会增加信号的长度，由卷积的运算过程可知，其中信号的头部与尾部是无用信息，因此需要把这部分无用信息去除掉，可以利用如下语句实现:

y_temp = conv(h,xdeal);                    

t = (N-1)/2;                          %N为滤波器的阶数

y = y_temp((t+1):(x_len+t));

3.8 去噪后音频的存盘

利用如下语句进行声音信号的读入：

[z,fs,nbits] = wavwrite(‘D:\...\...\result3.wav’);

  其中z为带存盘数据，fs为采样速率，nbits为数据z的位宽。

3.9 噪音区间的定位

（1）相关运算：

相关运算的过程相当于将处理过程中的一函数g经平移后，对另一函数f进行扫描；反应的是两函数在不同偏移量下的相似程度，相似程度的大小由它们的互相关函数值来衡量，在两函数相似程度最大的地方相关函数值最大。

相关系数的计算：

     X为模板信号的幅度，Y是待测信号的幅度。

  将上述相关系数计算公式编写成一个m文件，供主程序调用。

function coeff = correlation_ratio_calculation(X ,Y)

fenzi = sum(X .* Y) - (sum(X) * sum(Y)) / length(Y); 

fenmu = sqrt((sum(X .^2) - sum(X)^2 / length(Y)) * (sum(Y .^2) - sum(Y)^2 /     length(Y))); 

coeff = fenzi / fenmu; 

end

(2) 正弦、余弦模板函数的制作

N1 = 100;                  

y1 = zeros(1,N1);

y2 = zeros(1,N1);

f1 = 3500;                        %待检测噪音的频率

for i=1:1:N1

    y1(i) = sin(2*pi*i*f1/fs);    %频率为f1=3500的正弦信号模板

    y2(i) = cos(2*pi*i*f1/fs);    %频率为f1=3500的余弦信号模板

end

figure(5);

subplot(2,1,1);

plot(y1);

title('正弦信号模板');

subplot(2,1,2);

plot(y2);

title('余弦信号模板');

模板图形如下：

 

(3) 噪音的定位函数

%%%%利用相关系数定位噪音出现的时间段

ysin = zeros(1,x_len-N1);

ycos = zeros(1,x_len-N1);

rt   = zeros(1,x_len-N1);

for i=1:1:(x_len-N1)

    xt=xdeal(i:(N1+i-1));

    ysin(i) = correlation_ratio_calculation(xt,y1);

    ycos(i) = correlation_ratio_calculation(xt,y2);

    rt(i)   = abs(ysin(i))+abs(ycos(i));

end

   

figure(6);

plot(rt);

title('正弦余弦互补合成的相关系数');

figure(7);

subplot(2,1,1);

plot(ysin);

title('正弦模板相关系数');

subplot(2,1,2);

plot(ycos);

title('余弦模板相关系数');

噪声f1的定位：

有图像知，频率为f1 = 3500Hz的噪音出现的区间为[9900,196400].

噪音f2的定位：

有图像知，频率为f2 = 4600Hz的噪音出现的区间为[60000,196500].

噪音f3的定位：

有图像知，频率为f3 = 5500Hz的噪音出现的区间为[80000,196000].

4.实验结果与分析

用音乐播放器打开处理后的声音，可以听到噪音在一定程度上得到了抑制，能听出有用信息。但是与原始声音相比，效果还差很多。在同一幅图片中绘制夹杂噪音的信号与原始信号的频谱图：

对比两幅图可知，与原始信号相比噪音信号中除了3个正弦噪音外，还夹杂了许多其他频率的干扰，而所设计的滤波器只对3个正弦噪音做了处理，却没有处理其他干扰，这可能是滤波效果不太理想的一个重要原因。

5.总结与体会

通过这次课程设计，我不仅对matlab软件有了进一步的熟悉，更重要的是学会了如何对信号进行频域分析来找出信号中夹杂的噪音，以及如何设计滤波器来滤除噪音。还有如何利用相关算法来定位噪音在信号中出现的区间，体会到了学习《数字信号处理》重在应用与实践。

参考文献

1. 《数字信号处理》   程佩青著   清华大学出版社

2. 《MATLAB程序设计与应用基础教程》  张岳著  清华大学出版社

 

第二篇：数字信号处理课程设计三

课题一心电信号分析系统的设计

一、本课题的目的

本设计课题主要研究数字心电信号的初步分析方法及滤波器的应用。通过完成本课题的设计，拟主要达到以下几个目的：

（1）了解基于LabVIEW的虚拟仪器的特点和使用方法，熟悉采用LabVIEW进行仿真的方法。

（2）了解人体心电信号的时域特征和频谱特征。

（3）进一步了解数字信号的分析方法；

（4）通过应用具体的滤波器进一步加深对滤波器的理解。

（5）通过本课题的设计，培养学生运用所学知识分析和解决实际问题的能力。

二、课题任务

利用labVIEW设计一个基于虚拟仪器的简单的心电信号分析系统。对输入的原始心电信号，进行一定的数字信号处理，进行频谱分析。根据具体设计要求完成系统的程序编写、调试及功能测试。

（1）对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形。

（2）对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。

（3）根据心电信号的频域特征（自己查阅相关资料），设计相应的低通和带通滤波器。

（4）编程绘制实现信号处理前后的频谱，做频谱分析，得出相关结论。

（5）对系统进行综合测试，整理数据，撰写设计报告。

三、主要设备和软件

（1）PC机一台。

（2） LabVIEW软件一套，要求最低版本8.20。

四、设计内容、步骤和要求

必做部分：

1. 利用labVIEW读取MIT-BIH数据库提供的数字心电信号，并还原实际波形

美国麻省理工学院提供的MIT-BIH数据库是一个权威性的国际心电图检测标准库，近年来应用广泛，为我国的医学工程界所重视。MIT-BIH数据库共有48个病例，每个病例数据长30min，总计约有116000多个心拍，包含有正常心拍和各种异常心拍，内容丰富完整。

为了读取简单方便，采用其txt格式的数据文件作为我们的原心电信号数据。利用labVIEW提供的文件I/O函数，读取txt数据文件中的信号，并且还原实际波形。

2.对原始心电信号做线性插值处理

由于原始心电信号数据不是通过等间隔采样得到的，也就是说原始的心电数据并不是均匀的，而用Matlab中提供的数字滤波器处理数据时，要求数据是等间隔的。因此设计的系统首先应对原始心电信号做线性插值处理，使其变为等间隔的数字信号，否则直接处理后会出现偏差，根据心电信号的特点, 把时间分隔成0.001s。添加的幅值点采用一次线性插值。对二维数据进行插值,相连幅值间数据的插值根据时间进行,运算公式如下:

，，，，

其中是第i个数据时间点，A_i是与之对应的数据，N是两数据之间需要的插值数，是需要插值的两点数据差，，

时数组依次排列，即得到了插值后等间隔的新数据。

3.根据心电信号的频域特征，设计相应的低通和带通滤波器

一般正常人的心电信号频率在0.7～100HZ范围内，幅度为(胎儿)～5mV(成人)。人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴有干扰。采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz到0.25Hz频段的干扰，对于这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25,来滤除低频信号，对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99Hz。也可以直接应用带通滤波器设计。

根据以上说明，利用labVIEW中的信号处理函数设计相应滤波器，滤除数字信号中的干扰信号。

4.对处理前后的心电信号分别做频谱分析，分析结果

利用labVIEW对处理前后的心电信号编程显示其频谱，分析比对滤波处理前后的频谱，得出结论。

如果分析频谱，滤波效果不明显，则需变动滤波器参数指标，重新设计滤波器。通过频谱分析，多次试验确定最合适的滤波器。

5.系统界面设计

综合前面几步，设计出一个完整的系统，并且本着简洁的原则，设计友好的人机交互界面。

选作部分：

1.三种滤波器设计

分别设计Butterworth、Chebyshev、Inverse Chebyshev三种滤波器，并对滤波后的信号分别做频谱分析，比较几种滤波器的差别。

2.设计50HZ工频陷波器

由于电子设备采集到的信号经常会混有电源线干扰。电源线干扰是以50 Hz为中心的窄带噪声，带宽小于1Hz。设计相应滤波器滤除电源线干扰，并对处理后的信号做频谱分析。

五、课程设计报告要求

（1）设计报告书包括内容：课程设计题目，课程设计目的和意义，设计方案，详细设计步骤，设计结果（原理图等），测试和仿真结果（图形或数据）及其分析，其它有明确要求的设计内容，结论，参考文献等。

（2）提交课程设计报告时应同时提交相关设计和仿真分析材料（框图、程序、结果等）的电子版。

六、参考文献

[1] 陈锡辉，张银鸿编著.LabVIEW 8.20 程序设计从入门到精通[M].北京：清华大学出版社，2007.

[2] 丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，2001

[3] 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,8

[4] 谢嘉奎. 电子线路--线性部分(第四版). 高等教育出版社，2003,2

[5] 陈后金. 信号分析与处理实验. 高等教育出版社，2006，8

七、  附录——设计原理

附录：设计原理

1.心电信号的读取

txt格式的数据文件内容及格式如图1-1所示（以100.txt为例）：

图1.1 txt格式心电数据文件

其中文件的第一列为采样时间，第二列是在以MLII这种导联方式所得到的采样数据，第三列式以V5这种导联方式所得到的采样数据，全文件记录了约为10s的心电数据，3600个采样数据，每一行数据之间用Tab符分隔。

由于数据文件中后两列数据是对同一种心电信号进行不同的导联方式所得到的采样数据，所以可以采用任意其中的一种采样数据（比如选择MLII），摒弃另外一种，即可完成对此心电信号的分析。全部的心电文件记录时间约为10s，共计12个左右周期的心电信号。

根据txt格式的数据文件的特点，利用labvIEW提供的I/O文件函数，在本课题中，主要是围绕LabVIEW中的read from spreadsheet file读表单文件函数来设计心电信号的读取部分的VI，并利用XY Graph来对数据做图形化显示。让心电数据文件中的第一列时间数据作为x轴，对应的MLII方式的幅值作为y轴，以此得到绘制的原心电波形。

图1.2 读表单文件函数VI

图1.3读表单文件函数使用举例

实际设计心电信号数据文件时需要注意：

（1）数据文件的前两行为解释说明文字，不是真正的信号数据，读取信号程序要能够自动忽略前两行文字，只读取真正的数字信号数据（严禁自己手动删除原心电数据文件中的前两行数据，必须通过程序来实现忽略前两行文字的目的）。

（2）利用数组函数分别将文件的前两列分别读入一个一维数组。labvIEW默认的从文本文件中读取的数据都是字符串，因此在使用心电信号数据前需要将其转换为数值才可以。注意：第一列时间数据均为0:00.007这种格式，因此需要将字符串0:00.007先转化为字符串0.007，即去除字符串中冒号（：）以前的部分，然后再将其转为数值。

（3）最后利用已经转为数值的分别代表心电信号时间和幅值的两个一维数组，图形化还原原始心电信号波形，在此推荐利用labvIEW中XY Graph。

2.心电信号的线性插值处理

根据上文中提到的插值公式，以此为原理，设计labvIEW程序，对心电信号数据做线性插值处理。插值完以后的数据应该是时间均匀的、以0.001秒为间隔的。

此步骤主要是基于labvIEW中的数组操作函数来实现，建议一定首先熟悉并掌握labvIEW中的所有数组操作函数的作用和操作方法（比如array size函数、index array函数、insert into array函数等）。

其中一种插值方法的思路是：第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用for循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为0.001s，如果是则判断下一对相邻两个数据；如果间隔大于0.001s则在一个CASE结构里面做插值处理。

注意对时间数据做插值的同时一定不要忘记对幅值数据同样做插值处理，时间数据和幅值数据一定是相互对应的。

3.设计相应的数字滤波器

原心电信号里面是包含有噪声的，因此需要对数字心电信号做一定滤波处理。

LabVIEW提供的IIR滤波器类型有Butterworth、Chebyshev、Inverse Chebyshev、Elliptic和Bessel滤波器。它们都有各自的特点，用途也不尽相同。

LabVIEW还提供了高级IIR和FIR滤波器子面板。在高级面板中，滤波器的设计部分和执行部分是分开的。由于滤波器的设计很费时间，而滤波过程则很快。在含有循环结构的程序中，可以将滤波器的设计放在循环外，将设计好的滤波器参数传递到循环内，在循环内进行滤波，从而提高程序的运行效率。

Labview提供的滤波器函数面板面板位于Functions Palette的Signal Processing| Filters面板下，如图3.1所示。

图3.1 滤波器函数面板

选择合适的滤波器为心电信号设计一个低通和高通滤波器，或者带通滤波器。

4. 频谱分析

应该对线性插值后的心电信号和滤波处理后的心电信号做傅里叶变换，画出其频谱，比对前后差异，分析滤波器性能。

labvIEW中频域分析函数被划分为两个面板：Transforms面板实现的函数功能主要有傅立叶变换、Hilbert变换、小波变换、拉普拉斯变换等；Spectral Analysis面板包含的函数主要包括功率谱分析、联合时频分析等。

图4.1 transforms面板

 

图4.2 spectral analysis面板

5.低通滤波器和FFT举例

信号源由一个正弦信号与一个经过高通滤波的高频信号叠加而成。高通滤波器的截止频率为100Hz，即滤掉频率小于100 Hz的低频噪声分量。信号滤波器为Butterworth滤波器，截止频率设为30Hz，即滤掉频率大于30Hz的噪声分量。从图中可以清楚地看到滤波后的信号基本还原了正弦信号。

图4.3 低通滤波VI程序面板

图4.4 低通滤波VI前面板