信号与系统课程设计报告管道泄漏检测定位应用中的信号分析与处理Word格式.docx
《信号与系统课程设计报告管道泄漏检测定位应用中的信号分析与处理Word格式.docx》由会员分享，可在线阅读，更多相关《信号与系统课程设计报告管道泄漏检测定位应用中的信号分析与处理Word格式.docx（14页珍藏版）》请在冰点文库上搜索。
某学校有段自来水管道泄漏，该管道为铸铁管道，管径为250mm。
在管道两端接头处分别放置一个传感器及采集模块进行采集管道振动数据进行分析处理，得到漏点位置。
已知两个采集器的采样率均为10KHz，测试点之间的管道总长度L=100m，水声在管径为250mm铸铁管道中的传播速度v=1.250m/ms。
试图采用互相关方法对泄漏信号进行分析得到漏点位置，检测与定位原理如图1：
图1管道泄漏检测与定位原理
注：
L1,L2=漏点位置离传感器1、2的距离;
L=两传感器间的管道长度;
v=泄漏噪声在管道内的传播速度;
D=泄漏声信号到达两传感器的时间差.
泄漏点定位公式：
【设计要求】
1）对两端采集的数据分别进行时频域分析；
2）对两端采集的数据进行互相关运算，计算两路信号的相对延迟时间，并根据已给的两传感器间隔距离及水声沿管道传播速度计算泄漏点位置。
【设计步骤及提示】
通过MATLAB的函数load（）可以读入两个.txt格式的数据文件，并将该文件保存到指定的数组中。
例如下面的语句（更详细的命令介绍可以自己查阅MATLAB的帮助）中，将.txt读入后分别存放到矩阵y1和y2中。
y1=load（'
D:
\MATLAB\R2006a\work\S&
Sdesign\2005年12月2日15点39分set2daq1\data1.txt'
）;
y2=load（'
1.对泄漏信号进行时频域分析。
要求编写一个简单的函数文件SignalTF.m，它具有如下性质：
a、函数实现对信号进行快速傅里叶变换（FFT）；
b、绘制时域波形与频谱图。
要求图形中坐标标注齐全：
1）时域波形坐标标注为Time（s）—Amplitude（v），标题为“waveform”2）频谱图坐标标注为Frequency（Hz）—Magnitude（v/Hz）（即，横坐标标注为Frequency（Hz），纵坐标标注为Magnitude），标题为“Frequencyspectrum”。
找到泄漏信号的主要频谱成分所在的频带，验证为何检漏系统可以对泄漏信号采用10KHz的采样速率，并指出最低的采样频率是多少。
2.利用MATLAB函数库中函数设计butterworth高通滤波器（见Matlab的butter,filter函数），对采集的两路信号分别进行高通滤波，抑制工频噪声（典型值50Hz）及其他常见低频干扰噪声。
要求编写PreProcess.m函数，a）通过画图对滤波前后信号进行对比，观察差异，并用subplot命令分别将滤波前后的信号画在一个FIGURE图中（两个信号均需给出其滤波前后图形）；
b）将滤波后的两组数据分别写至data1p.txt和data2p.txt中。
3.编写TDE.m函数，用互相关方法检测两路信号的相对时间延迟，并以此计算漏点位置。
（提示：
用Matlab中xcorr函数实现两个信号互相关运算，根据峰值位置确定两路信号的相对时间间隔即数据点，并根据采样频率计算延迟时间，以此计算漏点位置。
）要求：
a）画出两路信号的互相关曲线图，找出最大点或峰值点及其横坐标，并标注齐全，峰值位置需对应正确的延迟时间（提示：
用Matlab中find函数查找最大值点及其位置）。
信号说明：
data1.txt：
传感器1采集的泄漏信号；
data2.txt：
传感器2采集的泄漏信号。
【MATLAB参考命令】
绘图命令：
plot、subplot等
求频谱、频率响应：
fft、freqz等
卷积、滤波：
conv、deconv、filter等
文件读写：
load、fopen,fprintf,fclose,xlsread,xlswrite等
三、正文
1.对泄露信号进行时频域分析
在MATLAB中，通过stem函数，可以画出离散时间信号在时域上的图像。
通过fft函数实现信号的快速傅里叶变换，使得信号从时域变换到频域上。
使用fftshift函数实现频率搬移，使频谱图关于原点左右对称。
这道题目的难点是确定频谱图的横坐标（即频率）的范围。
由采样定理，信号的最高频率不得超过采样频率的一半，否则采样后的信号将会出现失真。
由此可以确定频谱图的频率范围为-fs/2到fs/2（fs为采样频率）。
用函数linspace（-fs/2,fs/2,num）实现（num为采样信号的个数）
得到两个信号的时域图像如下：
得到信号的频谱图分别如下：
从频谱图可以看出，两个信号的最高频率大约为2200Hz，最低采样频率Fsmin=2*fm=4400Hz
程序代码如下：
functionY=SignalTF（X,N）%N为采样周期
n=N:
N:
length（X）*N;
%得到离散时间序列
fs=1/N;
num=length（X）;
Y=fft（X）;
Y1=fftshift（Y）;
figure,
stem（n,X,'
filled'
）%画出时域图像
title（'
waveform'
）
xlabel（'
Time（s）'
ylabel（'
Amlitude（v）'
mag=abs（Y1）;
%幅度
f=linspace（-fs/2,fs/2,num）;
%频率
plot（f,mag）%画频谱图
Frequecyspectrum'
Frequecy（Hz）'
Magnitude（v/Hz）'
2.
（1）利用MATLAB函数库中函数设计butterworth高通滤波器，对采集的两路信号分别进行高通滤波，抑制工频噪声（典型值50Hz）及其他常见低频干扰噪声。
MATLAB中的函数butter（N,Wn,’high’）实现butterworth高通滤波器，其中N表示滤波器的阶数，Wn表示截止频率。
工频噪声为50Hz，由于巴特沃兹滤波器是非理想滤波器，并且有其他频率的低频噪声，所以截止频率应高于50Hz，故这里截止频率取100Hz。
其中Wn应为0到1之间的数，Wn=截止频率/采样频率的一半。
通过filter这一维数字滤波器实现输出。
Y=filter（B,A,X），输入X为滤波前序列，Y为滤波结果序列，B/A提供滤波器的传递函数，B为分子，A为分母整个滤波过程是通过下面差分方程实现的：
a
（1）*y（n）=b
（1）*x（n）+b
（2）*x（n-1）+...+b（nb+1）*x（n-nb）-a
（2）*y（n-1）-...-a（na+1）*y（n-na）
[Y,Zf]=filter（B,A,X,Zi），输入X为滤波前序列，Y为滤波结果序列，B/A提供滤波器系数，B为分子，A为分母。
需要注意的是，这里的X和Y均为时域上的信号。
用subplot使滤波前后的两幅图画到同一张图上。
滤波前后的频谱图如下：
（上图为原始信号，下图为滤波后的信号）
从频谱图对比，可以看出，信号的低频分量已经被过滤掉了，滤波器设计成功。
具体代码如下：
functionY=PreProcess（X）%X为原始时域信号
[B,A]=butter（9,0.02,'
high'
%截止频率100Hz
F=fft（X）;
Y=filter（B,A,X）;
%其中X,Y均为时域信号
Z=fft（Y）;
Y1=fftshift（Z）;
fs=1e4;
%采样频率
subplot（2,1,1）%将滤波前后的信号画到一个图上
plot（f,abs（fftshift（F）））%滤波前的频谱图
xlabel（‘Frequency（Hz）’）
subplot（2,1,2）
plot（f,abs（Y1））%滤波后的频谱图
（2）将滤波后的两组数据分别写至data1p.txt和data2p.txt中。
MATLAB的文件操作跟C语言类似，用fopen打开一个文件，fprintf把数据写入文件，fclose关闭文件。
为了方便，这部分代码写在主程序下：
（主程序的完整代码见附录）
%以下分别将滤波后的两组数据写入data1p.txt和data2p.txt
fid1=fopen（'
H:
\学习\2012信号与系统-设计材料及要求\data1p.txt'
'
w'
fprintf（fid1,'
%12.8f\r\n'
F1）;
%记事本中，只能用\r\n才能换行
fclose（fid1）;
fid2=fopen（'
\学习\2012信号与系统-设计材料及要求\data2p.txt'
fprintf（fid2,'
F2）;
fclose（fid2）;
具体截图如下：
4.用互相关方法检测两路信号的相对时间延迟，并以此计算漏点位置
互相关函数表示的是两个时间序列之间和同一个时间序列在任意两个不同时刻的取值之间的相关程度，即互相关函数是描述随机信号x（t）,y（t）在任意两个不同时刻t1，t2的取值之间的相关程度。
对于连续信号公式表示为
R（τ）=（1/T）∫[f（t）g（t+τ）]dt，积分限为0至T。
对于离散信号公式表示为
R（n）=（1/N）∑[x（m）y（m+n）]其中m从0到N-1变化。
MATLAB中，用互相关函数xcorr，画出函数图象。
用max函数找出互相关函数的最大值所在的点。
设原信号的长度为N，采样频率为fs，最大值点为Nmax，则时间差D=（Nmax-N）/fs。
可以得出泄漏点距离第一个传感器的距离L1=（L-v*D）/2，到第二个泄漏点的距离为L2=（L+v*Dt）/2。
两个信号y1和y2的互相关函数的图象如下：
得
时间差D=0.0121s
泄漏点到第一个传感器的距离L1=42.4375m
泄漏点到第二个传感器的距离L2=57.5625m
functionY=TDE（x1,x2）
Y=xcorr（x1,x2）;
Lx=length（x1）;
Ly=length（Y）;
figure
plot（Y）,holdon
[Y_max,N]=max（Y）;
%求出最大值及对应的横坐标值
plot（N,Y_max,'
ro'
）%标出最大值点
Dn=N-Lx
D=Dn/fs%延迟时间
L=100;
v=1250;
L1=（L-v*D）/2
L2=（L+v*D）/2
四、结论
利用MATLAB可以对信号进行时域和频率域的分析，并且画出其时域图象和频谱图。
实际情况下找不到像理论课上所讲的理想滤波器，因此，在实际情况下使用butterworth高通滤波器可以很好地过滤掉低频噪声，并且阶数越高，越接近理想滤波器，滤波效果越好。
需要注意的是filter函数中，输入和输出信号都是在时域上的，否则会出现问题。
通过找出两个信号的互相关函数的最大值点，可以求出两个传感器接收信号的时间差，从而求出泄漏点到两个传感器之间的距离。
通过这次课程设计，我更加熟练地使用了MATLAB，对信号与系统的概念和物理意义有了更加深刻的理解，锻炼了我分析问题和解决问题的能力。
五、参考文献
[1].AlanV.Oppenheim,AlanS.Willsky,刘树棠译.信号与系统.西安交通大学出版社.2008.
[2].JohnR.Buck,MichaelM.Daniel,AndrewC.Singer,刘树棠译.信号与系统计算机练习——利用MATLAB.西安交通大学出版社.2006
六、附录
由于MATLAB不能像C语言那样，可以把多个函数写到同一个程序里。
故这里需要写一个主程序。
主程序代码如下：
%信号课程设计主程序
clear
y1=load（'
\学习\2012信号与系统-设计材料及要求\data1.txt'
y2=load（'
\学习\2012信号与系统-设计材料及要求\data2.txt'
f=1e4;
%采样率
N=1/f;
%采样周期
FT1=SignalTF（y1,N）;
F1=PreProcess（y1）;
FT2=SignalTF（y2,N）;
F2=PreProcess（y2）;
%记事本中，只能用\r\n才能换行
%第3问
TDE（y1,y2）;
升级会员 