数字信号课程设计空气主频率柱 洛阳理工

1、洛 阳 理 工 学 院课 程 设 计 报 告 课程名称 数字信号处理 设计题目 空气柱主频率模型测定 专 业 通信工程 班 级 B120509 学 号 B12050923 姓 名 徐 国 旗 完成日期 2015年6月27日 课 程 设 计 任 务 书设计题目： 空气柱主频率模型测定 _设计内容与要求： 题目3：空气柱主频率测定当我们向暖水瓶倒水时，暖水瓶中的空气柱会发生振动。随着水位的升高,空气柱越来越短，空气柱振动的频率越来越大，音调就越来越高，所以根据声音音调的高低就能知道水是否灌满。甚至我们可以量化这个过程，根据倒水声音音调的高低来反算水位的高低。笔者对倒水的过程进行了采集试验，选用标准

2、暖水瓶一个，用近似匀速的水流在43秒时间内将水瓶贮满，声音文件的STFT频谱图如下图所示：题目的基本要求如下：（1） 建立水位高度H和振动频率f的数学模型；（2） 自行设计一个倒水的过程，采集倒水过程的音频信号；（3） 用试验五的STFT程序对倒水声音进行分析；（4） 分离信号基本频率和水位高度的函数曲线；（5） 用函数曲线来验证数学模型的正确性。 课 程 设 计 评 语 成绩： 指导教师： 年 月 日洛 阳 理 工 学 院 课 程 设 计 报 告摘要数字信号处理课程是通信专业的一门重要专业基础课，是信息的数字化处理、存储和应用的基础。数字信号处理的主要研究对象是数字信号，且是采用运算的方法达

3、到处理的目的的，因此，其实现方法，基本上分成两种实现方法，即软件和硬件实现方法。软件实现方法指的是按照原理和算法，自己编写程序或者采用现成的程序在通用计算机上实现，硬件实现指的是按照具体的要求和算法，设计硬件结构图，用乘法器加法器延时器、控制器、存储器以及输入输出接口部件实现的一种方法。 语音信号的处理与滤波的设计主要是用Matlab作为工具平台，设计中涉及到声音的录制、播放、存储和读取，语音信号的抽样、频谱分析，滤波器的设计及语音信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来

4、达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。Matlab语言是一种广泛应用于工程计算及数值分析领域的新型高级语言，Matlab功能强大、简单易学、编程效率高，深受广大科技工作者的喜爱。特别是Matlab还具有信号分析工具箱，不需具备很强的编程能力，就可以很方便地进行语音信号分析、处理和设计。语音信号的处理与滤波是数字信号处理课程中常出现的课题，也是现代科学中值得深入究研的一个课题。 关键词：滤波器；信号处理；时频分析；数学建模目 录摘要1第1章 概述31.1 设计目的31.2 设计步骤31.2 MATLAB简介3第

5、2章 声音信号时域频域分析52.1时域频域波形显示源代码52.2时域频域波形效果图62.3时域频域结果分析6第3章 声音信号时频谱分析73.1时频谱显示源代码73.2时频谱显示效果图83.3时频谱显示结果分析8第4章 对声音信号进行滤波94.1 对原始信号进行切割94.1.1对原始信号进行切割并显示源代码94.2 对切割后信号分段进行滤波104.2.1对切割后信号进行滤波并显示源代码104.2.2滤波器特性及滤波后显示效果图124.2.3切割并滤波后显示效果分析124.3对滤波后的信号进行合成134.3.1 对滤波后声音信号进行合成源代码134.3.2.合成后显示效果图13第5章 从合成信号中

6、分离频率与时间的曲线145.1分离频率与时间的关系程序代码145.2分离后的频率与时间的关系图145.3对分离后的频率与时间的关系分析14第6章 建立频率与水位高度的数学模型156.1建立频率与时间的数学模型156.2频率与时间的数学模型结果图156.3建立频率与水位高度的函数模型16心得体会17参考文献18第1章 概述1.1 设计目的（1）学会MATLAB的使用，掌握MATLAB的程序设计方法；（2）掌握数字信号处理的基本概念、基本理论和基本方法；（3）学会用MATLAB对音频信号进行分析和处理；（4）掌握MATLAB设计滤波器；（5）学会用MATLAB的工具箱进行数学建模，解决一些实际问题

7、。1.2 设计步骤（1）对倒水语音信号进行采集；（2）用MATLAB对倒水语音信号进行时域频域分析；（3）用MATLAB对倒水语音信号进行时频联合分析；（4）对原始声音信号进行切割并分段滤波；（5）对滤波后的信号进行合成；（6）分离信号函数曲线；（7）建立数学模型。1.2 MATLAB简介MATLAB是矩阵实验室(Matrix Laboratory)之意。除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。MATLAB的基本数据单位是矩阵，它的指令表达式与数学，工程中常用的形式十分相似，故用MATLAB来解算问题要比用C,FORTRAN等语言完相

8、同的事情简捷得多，工具包又可以分为功能性工具包和学科工具包。功能工具包用来扩充MATLAB的符号计算，可视化建模仿真，文字处理及实时控制等功能。学科工具包是专业性比较强的工具包，控制工具包，信号处理工具包，通信工具包等都属于此类。MATLAB具有许多的优点比如：语言简洁紧凑，使用方便灵活，库函数极其丰富；MATLAB既具有结构化的控制语句(如for循环，while循环，break语句和if语句)，又有面向对象编程的特性;程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行，等等优点。MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建

9、模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用MATLAB 函数集）扩展了MATLAB 环境，以解决这些应用领域内特定类型的问题。其启动界面如图1-1所示。图1-1 MATLAB启动界面第2章 声音信号时域频域分析2.1时域频域波形显示源代码clearcloseclcy,fs=wavread('daoshui.wav');% 显示声音信号时域波形t = linspace(0,50,length(y);subplot(211);plot(t,y);%添加网格，标题和坐标轴title('声音信号的时域波形');xlabel('时间(s)&

10、#39;);ylabel('幅度');% 显示声音信号频域波形fsig=abs(fftshift(fft(y); ww = linspace(-fs/2, fs/2, length(fsig); subplot(212);plot(ww,20*log10(fsig);%添加网格，标题和坐标轴grid on;title('声音信号的频域波形');xlabel('频率（HZ）');ylabel('幅度');2.2时域频域波形效果图图2-1时域频域波形2.3时域频域结果分析从时域频域谱中只能看出含有各个频率分量，但并不能看出具体有什么各

11、频率成分，所以需要时频联合分析。第3章 声音信号时频谱分析3.1时频谱显示源代码Clear close clcy,fs = wavread('daoshui.wav');%持续时间T = length(y)/8000;MF = 4000;%对信号进行时频分析 I t f = ShowSound(y,8000,1.6);%构造时间频率轴和时频谱缓冲区t = linspace(0,T,size(I,2);f = linspace(0,MF,size(I,1);II = zeros(size(I);II(end:-1:1,:) = I(:,:);mesh(t,f,log2(II);t

12、itle('原始信号时频谱');xlabel('时间(s)');ylabel('频谱(Hz)')zlabel('幅度(dB)');view(0,90); 程序中调用的ShowSound函数是老师给定的，其代码如下：function Y,T,F = ShowSound(s,fs,lamda)%用矩形窗的FFT将采样率为fs的一段声音s进行时频分析%Y为输出的时频图样，T是时间轴，F频率轴；%由于显示设备的限制，时间轴只能显示750个点，频率轴能显示500个点%lamda是矩形窗的交叠系数，即交叠的部分占总窗体的百分比%作者：王新新

13、时间：2005.8.2%统一为行向量r,c = size(s);if(r > c) s = s'end;LS=length(s);%加750次窗，每次增加的点数SubLen = floor(LS / 750);%将序列化为750的整倍数s = s(1 : 750 * SubLen);LS = length(s);%由参数lamda计算出真实窗体的长度TrueSubLen = floor(SubLen./ (1 - lamda./ 2);%将序列扩展以便于连续加窗offset = TrueSubLen - SubLen;s = zeros(1,offset),s;%将加窗后的全部序

14、列暂存于tmp，用for循环生成tmptmp = zeros(TrueSubLen,750);for i = 2 : 750 index = floor(i - 1) * SubLen - (TrueSubLen * (lamda./ 2)+offset; tmp(:,i) = (s(index : (index+TrueSubLen-1)'end;%计算FFT并将颠倒的谱倒过来，只取频谱的正半边Y = fft(tmp , 1000);Y = Y(500:end , :);Y = (real(Y). 2 + imag(Y). 2).0.5;%计算T和FT = 0 : 1 / fs :

15、(length(s) - 1)./ fs;F = 0 : 500;3.2时频谱显示效果图图3-1 时频谱显示效果图3.3时频谱显示结果分析由图3-1时频谱可看出各频率分量是随时间变化的，而且后一时刻的频率可能会和前面某一时刻的频率大小相等，所以要想通过滤波取主频率柱，必须先对原始信号进行切割，滤波后再合成。第4章 对声音信号进行滤波4.1 对原始信号进行切割4.1.1对原始信号进行切割并显示源代码clear all close all clcy,fs=wavread('daoshui.wav');MF = 4000;I t f = ShowSound(y,8000,1.6);%

16、切割第1段0-15sT11=1/fs;T12 = 15;y1=y(T11*fs:T12*fs); %切割第2段15-25sT21=15+1/fs;T22=25;y2=y(T21*fs:T22*fs); %切割第3段25-30sT31=25+1/fs;T32=30;y3=y(T31*fs:T32*fs);%切割第4段30-35sT41=30+1/fs;T42=35;y4=y(T41*fs:T42*fs);%切割第5段35-44sT51=35+1/fs;T52=44;y5=y(T51*fs:T52*fs);4.2 对切割后信号分段进行滤波4.2.1对切割后信号进行滤波并显示源代码滤波时采用巴特沃夫

17、滤波器，其滤波器特性如下图4-1所示。%滤波并显示程序% 第1段切割并滤波后B,A = butter(4,600/(fs/2),1150/(fs/2), 'bandpass');res1=filter(B,A,y1);I t f = ShowSound(res1,8000,1.6);t = linspace(T11,T12,size(I,2);f = linspace(0,4000,size(I,1);II = zeros(size(I);II(end:-1:1,:) = I(:,:);subplot(231);mesh(t,f,log2(II);axis tighttitle

18、('第1段切割并滤波后');xlabel('时间(s)');ylabel('频率');view(0,90);% 第2段切割并滤波后B,A = butter(3,900/(fs/2),1800/(fs/2), 'bandpass');res2=filter(B,A,y2);I t f = ShowSound(res2,8000,1.6);t = linspace(T21,T22,size(I,2);f = linspace(0,4000,size(I,1);II = zeros(size(I);II(end:-1:1,:) = I(

19、:,:);subplot(232);mesh(t,f,log2(II);axis tighttitle('第2段切割并滤波后');xlabel('时间(s)');ylabel('频率');view(0,90);% 第3段切割并滤波后B,A = butter(3,1600/(fs/2),2200/(fs/2), 'bandpass');res3=filter(B,A,y3);I t f = ShowSound(res3,8000,1.6);t = linspace(T31,T32,size(I,2);f = linspace(0,4

20、000,size(I,1);II = zeros(size(I);II(end:-1:1,:) = I(:,:);subplot(233);mesh(t,f,log2(II);axis tighttitle('第3段切割并滤波后');xlabel('时间(s)');ylabel('频率');view(0,90);% 第4段切割并滤波后B,A = butter(3,2000/(fs/2),2600/(fs/2), 'bandpass');res4=filter(B,A,y4);I t f = ShowSound(res4,8000,

21、1.6);t = linspace(T41,T42,size(I,2);f = linspace(0,4000,size(I,1);II = zeros(size(I);II(end:-1:1,:) = I(:,:);subplot(234);mesh(t,f,log2(II);axis tighttitle('第4段切割并滤波后');xlabel('时间(s)');ylabel('频率');view(0,90);% 第5段切割并滤波后B,A = butter(3,2400/(fs/2),3200/(fs/2), 'bandpass

22、9;);res5=filter(B,A,y5);I t f = ShowSound(res5,8000,1.6);t = linspace(T51,T52,size(I,2);f = linspace(0,4000,size(I,1);II = zeros(size(I);II(end:-1:1,:) = I(:,:);subplot(235);mesh(t,f,log2(II);axis tighttitle('第5段切割并滤波后');xlabel('时间(s)');ylabel('频率');view(0,90);4.2.2滤波器特性及滤波后显

23、示效果图图4-1滤波器特性图图4-2切割并滤波后显示效果图4.2.3切割并滤波后显示效果分析由图可以看出，前几段滤波效果比较好，而最后两段滤波效果不是太明显，原因是在即将倒满水时频率分量较杂，所需主频率幅度较小。但可从前面的已知数据中通过数学手段预测出来。4.3对滤波后的信号进行合成4.3.1 对滤波后声音信号进行合成源代码%合成程序%yres(T11*fs:T12*fs)=res1;yres(T21*fs:T22*fs)=res2;yres(T31*fs:T32*fs)=res3;yres(T41*fs:T42*fs)=res4;yres(T51*fs:T52*fs)=res5;% 显示合成

24、效果I t f = ShowSound(yres,8000,1.6);t = linspace(T11,T52,size(I,2);f = linspace(0,MF,size(I,1);II = zeros(size(I);II(end:-1:1,:) = I(:,:);figuremesh(t,f,log2(II);axis tighttitle('合成效果图');xlabel('时间(s)');ylabel('频率');view(0,90);4.3.2.合成后显示效果图图4-3 滤波合成后显示效果第5章 从合成信号中分离频率与时间的曲线5.

25、1分离频率与时间的关系程序代码%根据幅度峰值找出频率F与时间t的关系%F = zeros(1,750);for T=1:750 m,F(T)=max(II(:,T);endt = linspace(0,43,750);F=(F/500)*4000; %显示范围从0-500变为0-4000figuretitle('频率F与时间t的关系图');xlabel('时间');ylabel('频率');plot(t,F);5.2分离后的频率与时间的关系图图5-1 频率与时间的函数曲线图5.3对分离后的频率与时间的关系分析通过图像可以看出，频率与时间的关系前面

26、效果还是比较好的，只是最后一小段不太好，但这并不影响进行数学建模。第6章 建立频率与水位高度的数学模型6.1建立频率与时间的数学模型%数学建模%yy=F(16:700);xx=t(16:700);figureplot(xx,yy,'r.');%先描点，看看大致是什么图，这里看到像指数函数曲线。hold on ;p=polyfit(xx,log(yy),3);xx=t(16:750);X=polyval(p,xx);Y=exp(X);plot(xx,Y,'.')title('频率和时间模型图');xlabel('时间');ylabe

27、l('频率');根据运行结果，p =-0.000018 0. 001183 0.015755 6.6419125可得到式6-1的函数模型。6-1X为时间，Y为频率。P为系数。6.2频率与时间的数学模型结果图图6-1 频率与时间的函数模型结果图6.3建立频率与水位高度的函数模型由于水是用近似匀速的水流在43秒时间内将水瓶贮满，而经测量，暖水壶高度近似于35cm。这里忽略暖水瓶最低处和最高处的弯曲，把其等效成圆柱体，所以可以用以下程序得出频率与水位高度的函数模型：figureh = linspace(0,35,length(xx);plot(h,Y,'.')titl

28、e('频率和水位高度模型图');xlabel('水位高度h(cm)');ylabel('频率F(HZ)');图6-2 频率F与水位高度h的函数模型图可以看出当水位不断升高时频率不断增大。函数模型如式6-2所示：6-2至此为止，已建立了频率F与水位高度h的数学模型。然后就可以根据此模型用单片机等制作通过测量频率计算暖水瓶水位高度的仪器，并加入其它模型，使其变成多用途测量仪，相比其它需要用直接接触的传感器来测量的仪器，有其独特的优势。心得体会这次课程设计令我印象非常深刻，通过这次课程设计，我不仅对自己的理论知识学习的更加扎实，而且锻炼了我的动手能力、设计创新的能力。这次课题布置的比较早，在布置之初我就开