信號(hào)與系統(tǒng)課程設(shè)計(jì)報(bào)告-管道泄漏檢測(cè)定位應(yīng)用中的信號(hào)分析與處理

1、<p>　　《信號(hào)與系統(tǒng)》課程設(shè)計(jì)報(bào)告</p><p>　　管道泄漏檢測(cè)定位應(yīng)用中的信號(hào)分析與處理</p><p>　　學(xué)院： </p><p>　　班級(jí)： </p><p>　　姓名： </p><p>　　學(xué)號(hào)：

2、 </p><p>　　指導(dǎo)老師： </p><p><b>　　本小組其他成員： </b></p><p><b>　　摘要</b></p><p>　　管道泄漏時(shí)，泄漏點(diǎn)將會(huì)產(chǎn)生振動(dòng)，出現(xiàn)聲音信號(hào)。</p><p>　　通過對(duì)兩個(gè)采樣器采集

3、到的振動(dòng)信號(hào)進(jìn)行時(shí)域和頻域的分析，使用巴特沃茲高通濾波器，濾掉低頻（50Hz）噪聲信號(hào)。采用互相關(guān)函數(shù)找到兩個(gè)采樣器接收信號(hào)的時(shí)間差，根據(jù)聲音在管道中的速度，求出管道泄漏點(diǎn)的位置。利用MATLAB實(shí)現(xiàn)以上功能。</p><p>　　關(guān)鍵字：時(shí)頻域分析、巴特沃茲高通濾波器、互相關(guān)函數(shù)</p><p><b>　　二、設(shè)計(jì)原理及要求</b></p><

4、p>　　【管道泄漏檢測(cè)原理】</p><p>　　已知：某學(xué)校有段自來水管道泄漏，該管道為鑄鐵管道，管徑為250mm。在管道兩端接頭處分別放置一個(gè)傳感器及采集模塊進(jìn)行采集管道振動(dòng)數(shù)據(jù)進(jìn)行分析處理，得到漏點(diǎn)位置。已知兩個(gè)采集器的采樣率均為10KHz，測(cè)試點(diǎn)之間的管道總長(zhǎng)度L=100m，水聲在管徑為250mm鑄鐵管道中的傳播速度v=1.250m/ms。試圖采用互相關(guān)方法對(duì)泄漏信號(hào)進(jìn)行分析得到漏點(diǎn)位置，檢測(cè)與定位

5、原理如圖1：</p><p>　　圖1 管道泄漏檢測(cè)與定位原理</p><p><b>　　注：</b></p><p>　　L1, L2 = 漏點(diǎn)位置離傳感器1、2的距離; </p><p>　　L = 兩傳感器間的管道長(zhǎng)度; </p><p>　　v = 泄漏噪聲在管道內(nèi)的傳播速度; &l

6、t;/p><p>　　D = 泄漏聲信號(hào)到達(dá)兩傳感器的時(shí)間差. </p><p><b>　　泄漏點(diǎn)定位公式：</b></p><p><b>　　【設(shè)計(jì)要求】</b></p><p>　　對(duì)兩端采集的數(shù)據(jù)分別進(jìn)行時(shí)頻域分析；</p><p>　　對(duì)兩端采集的數(shù)據(jù)進(jìn)行互相關(guān)運(yùn)算，

7、計(jì)算兩路信號(hào)的相對(duì)延遲時(shí)間，并根據(jù)已給的兩傳感器間隔距離及水聲沿管道傳播速度計(jì)算泄漏點(diǎn)位置。</p><p><b>　　【設(shè)計(jì)步驟及提示】</b></p><p>　　通過MATLAB的函數(shù)load ()可以讀入兩個(gè).txt格式的數(shù)據(jù)文件，并將該文件保存到指定的數(shù)組中。例如下面的語句（更詳細(xì)的命令介紹可以自己查閱MATLAB的幫助）中，將.txt讀入后分別存放到矩陣

8、y1和y2中。</p><p>　　y1 = load(' D:\MATLAB\R2006a\work\S&Sdesign\2005年12月2日15點(diǎn)39分set2daq1\data1.txt'); </p><p>　　y2 = load(' D:\MATLAB\R2006a\work\S&Sdesign\2005年12月2日15點(diǎn)39分set2

9、daq1\data1.txt'); </p><p>　　對(duì)泄漏信號(hào)進(jìn)行時(shí)頻域分析。要求編寫一個(gè)簡(jiǎn)單的函數(shù)文件SignalTF.m，它具有如下性質(zhì)：a、函數(shù)實(shí)現(xiàn)對(duì)信號(hào)進(jìn)行快速傅里葉變換(FFT)；b、繪制時(shí)域波形與頻譜圖。要求圖形中坐標(biāo)標(biāo)注齊全：1）時(shí)域波形坐標(biāo)標(biāo)注為Time（s）—Amplitude（v），標(biāo)題為“waveform”2）頻譜圖坐標(biāo)標(biāo)注為Frequency（Hz）—Magnitude（v

10、/Hz）（即，橫坐標(biāo)標(biāo)注為Frequency（Hz），縱坐標(biāo)標(biāo)注為Magnitude），標(biāo)題為“Frequency spectrum”。找到泄漏信號(hào)的主要頻譜成分所在的頻帶，驗(yàn)證為何檢漏系統(tǒng)可以對(duì)泄漏信號(hào)采用10KHz的采樣速率，并指出最低的采樣頻率是多少。</p><p>　　利用MATLAB函數(shù)庫中函數(shù)設(shè)計(jì)butterworth高通濾波器（見Matlab 的butter, filter函數(shù)），對(duì)采集的兩路信號(hào)

11、分別進(jìn)行高通濾波，抑制工頻噪聲（典型值50Hz）及其他常見低頻干擾噪聲。要求編寫PreProcess.m函數(shù)，a) 通過畫圖對(duì)濾波前后信號(hào)進(jìn)行對(duì)比，觀察差異，并用subplot命令分別將濾波前后的信號(hào)畫在一個(gè)FIGURE圖中（兩個(gè)信號(hào)均需給出其濾波前后圖形）；b) 將濾波后的兩組數(shù)據(jù)分別寫至data1p.txt和data2p.txt中。</p><p>　　編寫TDE.m函數(shù)，用互相關(guān)方法檢測(cè)兩路信號(hào)的相對(duì)時(shí)間延

12、遲，并以此計(jì)算漏點(diǎn)位置。（提示：用Matlab中xcorr函數(shù)實(shí)現(xiàn)兩個(gè)信號(hào)互相關(guān)運(yùn)算，根據(jù)峰值位置確定兩路信號(hào)的相對(duì)時(shí)間間隔即數(shù)據(jù)點(diǎn)，并根據(jù)采樣頻率計(jì)算延遲時(shí)間，以此計(jì)算漏點(diǎn)位置。）要求：a) 畫出兩路信號(hào)的互相關(guān)曲線圖，找出最大點(diǎn)或峰值點(diǎn)及其橫坐標(biāo)，并標(biāo)注齊全，峰值位置需對(duì)應(yīng)正確的延遲時(shí)間（提示：用Matlab中find函數(shù)查找最大值點(diǎn)及其位置）。</p><p><b>　　信號(hào)說明：</b&

13、gt;</p><p>　　data1.txt：傳感器1采集的泄漏信號(hào)；data2.txt：傳感器2采集的泄漏信號(hào)。</p><p>　　【MATLAB參考命令】</p><p>　　繪圖命令：plot、subplot等</p><p>　　求頻譜、頻率響應(yīng)：fft、freqz等</p><p>　　卷積、濾波：con

14、v、deconv、filter等</p><p>　　文件讀寫：load、fopen, fprintf, fclose, xlsread, xlswrite等</p><p><b>　　三、正文</b></p><p>　　1.對(duì)泄露信號(hào)進(jìn)行時(shí)頻域分析</p><p>　　在MATLAB中，通過stem函數(shù)，可以畫出離

15、散時(shí)間信號(hào)在時(shí)域上的圖像。通過fft函數(shù)實(shí)現(xiàn)信號(hào)的快速傅里葉變換，使得信號(hào)從時(shí)域變換到頻域上。使用fftshift函數(shù)實(shí)現(xiàn)頻率搬移，使頻譜圖關(guān)于原點(diǎn)左右對(duì)稱。這道題目的難點(diǎn)是確定頻譜圖的橫坐標(biāo)（即頻率）的范圍。由采樣定理，信號(hào)的最高頻率不得超過采樣頻率的一半，否則采樣后的信號(hào)將會(huì)出現(xiàn)失真。由此可以確定頻譜圖的頻率范圍為-fs/2 到 fs/2（fs為采樣頻率）。用函數(shù)linspace(-fs/2,fs/2,num)實(shí)現(xiàn)（num為采樣信號(hào)

16、的個(gè)數(shù)） </p><p>　　得到兩個(gè)信號(hào)的時(shí)域圖像如下：</p><p>　　得到信號(hào)的頻譜圖分別如下：</p><p>　　從頻譜圖可以看出，兩個(gè)信號(hào)的最高頻率大約為2200Hz，最低采樣頻率Fsmin=2*fm=4400Hz</p><p><b>　　程序代碼如下：</b></p><p&

17、gt;　　function Y=SignalTF(X,N) %N為采樣周期</p><p>　　n=N:N:length(X)*N; %得到離散時(shí)間序列</p><p><b>　　fs=1/N;</b></p><p>　　num=length(X);</p><p><b>　　Y=fft(X);&

18、lt;/b></p><p>　　Y1=fftshift(Y);</p><p><b>　　figure,</b></p><p>　　stem(n,X,'filled') %畫出時(shí)域圖像</p><p>　　title('waveform')</p><

19、;p>　　xlabel('Time(s)')</p><p>　　ylabel('Amlitude(v)')</p><p>　　mag=abs(Y1); %幅度</p><p><b>　　figure,</b></p><p>　　f=linspace(-fs/2,fs/2,nu

20、m); %頻率</p><p>　　plot(f,mag) %畫頻譜圖</p><p>　　title('Frequecy spectrum')</p><p>　　xlabel('Frequecy(Hz)')</p><p>　　ylabel('Magnitude(v/Hz)')<

21、/p><p>　　2.（1）利用MATLAB函數(shù)庫中函數(shù)設(shè)計(jì)butterworth高通濾波器，對(duì)采集的兩路信號(hào)分別進(jìn)行高通濾波，抑制工頻噪聲（典型值50Hz）及其他常見低頻干擾噪聲。</p><p>　　MATLAB中的函數(shù)butter(N,Wn,’high’)實(shí)現(xiàn)butterworth高通濾波器，其中N表示濾波器的階數(shù)，Wn表示截止頻率。工頻噪聲為50Hz，由于巴特沃茲濾波器是非理想濾波器，

22、并且有其他頻率的低頻噪聲，所以截止頻率應(yīng)高于50Hz，故這里截止頻率取100Hz。其中Wn應(yīng)為0到1之間的數(shù)，Wn=截止頻率/采樣頻率的一半。通過filter這一維數(shù)字濾波器實(shí)現(xiàn)輸出。Y = filter(B,A,X) ，輸入X為濾波前序列，Y為濾波結(jié)果序列，B/A 提供濾波器的傳遞函數(shù)，B為分子， A為分母 整個(gè)濾波過程是通過下面差分方程實(shí)現(xiàn)的： a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... +

23、 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為濾波結(jié)果序列，B/A 提供濾波器系數(shù)，B為分子， A為分母。需要注意的是，這里的X和Y均為時(shí)域上的信號(hào)。用subplot使濾波前后的兩幅圖畫到同一張圖上。</p><p>　　濾波前后的頻譜圖如下：(上圖為原始信號(hào)，下圖為濾波

24、后的信號(hào))</p><p>　　從頻譜圖對(duì)比，可以看出，信號(hào)的低頻分量已經(jīng)被過濾掉了，濾波器設(shè)計(jì)成功。</p><p><b>　　具體代碼如下：</b></p><p>　　function Y=PreProcess(X) %X為原始時(shí)域信號(hào)</p><p>　　num=length(X);</p>

25、<p>　　[B,A]=butter(9,0.02,'high');%截止頻率100Hz</p><p><b>　　F=fft(X);</b></p><p>　　Y=filter(B,A,X); %其中X,Y均為時(shí)域信號(hào)</p><p><b>　　Z=fft(Y);</b></

26、p><p>　　Y1=fftshift(Z);</p><p>　　fs=1e4; %采樣頻率</p><p><b>　　figure,</b></p><p>　　f=linspace(-fs/2,fs/2,num);</p><p>　　subplot(2,1,1) %將濾波前后的信號(hào)畫到一個(gè)

27、圖上</p><p>　　plot(f,abs(fftshift(F))) %濾波前的頻譜圖</p><p>　　xlabel(‘Frequency(Hz)’)</p><p>　　subplot(2,1,2)</p><p>　　plot(f,abs(Y1)) %濾波后的頻譜圖</p><p>　　xlabel

28、(‘Frequency(Hz)’)</p><p>　　（2）將濾波后的兩組數(shù)據(jù)分別寫至data1p.txt和data2p.txt中。</p><p>　　MATLAB的文件操作跟C語言類似，用fopen打開一個(gè)文件，fprintf把數(shù)據(jù)寫入文件，fclose關(guān)閉文件。為了方便，這部分代碼寫在主程序下：（主程序的完整代碼見附錄）</p><p>　　%以下分別將濾波

29、后的兩組數(shù)據(jù)寫入data1p.txt 和data2p.txt</p><p>　　fid1=fopen('H:\學(xué)習(xí)\2012信號(hào)與系統(tǒng)-設(shè)計(jì)材料及要求\data1p.txt','w');</p><p>　　fprintf(fid1,'%12.8f\r\n',F1); %記事本中，只能用\r\n才能換行</p><p&

30、gt;　　fclose(fid1);</p><p>　　fid2=fopen('H:\學(xué)習(xí)\2012信號(hào)與系統(tǒng)-設(shè)計(jì)材料及要求\data2p.txt','w');</p><p>　　fprintf(fid2,'%12.8f\r\n',F2);</p><p>　　fclose(fid2);</p>&

31、lt;p><b>　　具體截圖如下：</b></p><p>　　用互相關(guān)方法檢測(cè)兩路信號(hào)的相對(duì)時(shí)間延遲，并以此計(jì)算漏點(diǎn)位置</p><p>　　互相關(guān)函數(shù)表示的是兩個(gè)時(shí)間序列之間和同一個(gè)時(shí)間序列在任意兩個(gè)不同時(shí)刻的取值之間的相關(guān)程度，即互相關(guān)函數(shù)是描述隨機(jī)信號(hào)x(t),y(t)在任意兩個(gè)不同時(shí)刻t1，t2的取值之間的相關(guān)程度。對(duì)于連續(xù)信號(hào)公式表示為</p

32、><p>　　R(τ)=（1/T）∫[f(t)g(t+τ)]dt，積分限為0至T。對(duì)于離散信號(hào)公式表示為</p><p>　　R(n)=（1/N）∑[x(m)y(m+n)]其中m從0到N-1變化。</p><p>　　MATLAB中，用互相關(guān)函數(shù)xcorr，畫出函數(shù)圖象。用max函數(shù)找出互相關(guān)函數(shù)的最大值所在的點(diǎn)。設(shè)原信號(hào)的長(zhǎng)度為N，采樣頻率為fs，最大值點(diǎn)為Nmax，

33、則時(shí)間差D=(Nmax-N)/fs?？梢缘贸鲂孤c(diǎn)距離第一個(gè)傳感器的距離L1=（L-v*D）/2，到第二個(gè)泄漏點(diǎn)的距離為L(zhǎng)2=(L+v*Dt)/2。</p><p>　　兩個(gè)信號(hào)y1和y2的互相關(guān)函數(shù)的圖象如下：</p><p><b>　　得</b></p><p>　　時(shí)間差 D=0.0121s</p><p>　　

34、泄漏點(diǎn)到第一個(gè)傳感器的距離 L1= 42.4375m</p><p>　　泄漏點(diǎn)到第二個(gè)傳感器的距離 L2= 57.5625m</p><p><b>　　具體代碼如下：</b></p><p>　　function Y=TDE(x1,x2)</p><p>　　Y=xcorr(x1,x2);</p>&l

35、t;p>　　Lx=length(x1);</p><p>　　Ly=length(Y);</p><p><b>　　fs=1e4;</b></p><p><b>　　figure</b></p><p>　　plot(Y), hold on</p><p>　　[Y

36、_max,N]=max(Y); %求出最大值及對(duì)應(yīng)的橫坐標(biāo)值</p><p>　　plot(N,Y_max,'ro') %標(biāo)出最大值點(diǎn)</p><p><b>　　Dn=N-Lx</b></p><p>　　D=Dn/fs %延遲時(shí)間</p><p><b>　　L=100;<

37、/b></p><p><b>　　v=1250;</b></p><p>　　L1=(L-v*D)/2</p><p>　　L2=(L+v*D)/2</p><p><b>　　四、結(jié)論</b></p><p>　　利用MATLAB可以對(duì)信號(hào)進(jìn)行時(shí)域和頻率域的分析，并

38、且畫出其時(shí)域圖象和頻譜圖。實(shí)際情況下找不到像理論課上所講的理想濾波器，因此，在實(shí)際情況下使用butterworth高通濾波器可以很好地過濾掉低頻噪聲，并且階數(shù)越高，越接近理想濾波器，濾波效果越好。需要注意的是filter函數(shù)中，輸入和輸出信號(hào)都是在時(shí)域上的，否則會(huì)出現(xiàn)問題。通過找出兩個(gè)信號(hào)的互相關(guān)函數(shù)的最大值點(diǎn)，可以求出兩個(gè)傳感器接收信號(hào)的時(shí)間差，從而求出泄漏點(diǎn)到兩個(gè)傳感器之間的距離。</p><p>　　通過這

39、次課程設(shè)計(jì)，我更加熟練地使用了MATLAB，對(duì)信號(hào)與系統(tǒng)的概念和物理意義有了更加深刻的理解，鍛煉了我分析問題和解決問題的能力。</p><p><b>　　五、參考文獻(xiàn)</b></p><p>　　[1].Alan V. Oppenheim, Alan S. Willsky, 劉樹棠譯.信號(hào)與系統(tǒng).西安交通大學(xué)出版社.2008.</p><p>

40、;　　[2].John R. Buck, Michael M. Daniel, Andrew C. Singer, 劉樹棠譯.信號(hào)與系統(tǒng)計(jì)算機(jī)練習(xí)——利用MATLAB.西安交通大學(xué)出版社.2006</p><p><b>　　六、附錄</b></p><p>　　由于MATLAB不能像C語言那樣，可以把多個(gè)函數(shù)寫到同一個(gè)程序里。故這里需要寫一個(gè)主程序。主程序代碼如下：

41、</p><p>　　%信號(hào)課程設(shè)計(jì)主程序</p><p><b>　　clear</b></p><p>　　y1=load('H:\學(xué)習(xí)\2012信號(hào)與系統(tǒng)-設(shè)計(jì)材料及要求\data1.txt');</p><p>　　y2=load('H:\學(xué)習(xí)\2012信號(hào)與系統(tǒng)-設(shè)計(jì)材料及要求\data

42、2.txt');</p><p>　　f=1e4; %采樣率</p><p>　　N=1/f; %采樣周期</p><p>　　FT1=SignalTF(y1,N);</p><p>　　F1=PreProcess(y1);</p><p>　　FT2=SignalTF(y2,N);</p>

43、<p>　　F2=PreProcess(y2);</p><p>　　%以下分別將濾波后的兩組數(shù)據(jù)寫入data1p.txt 和data2p.txt</p><p>　　fid1=fopen('H:\學(xué)習(xí)\2012信號(hào)與系統(tǒng)-設(shè)計(jì)材料及要求\data1p.txt','w');</p><p>　　fprintf(fid1,&#

44、39;%12.8f\r\n',F1); %記事本中，只能用\r\n才能換行</p><p>　　fclose(fid1);</p><p>　　fid2=fopen('H:\學(xué)習(xí)\2012信號(hào)與系統(tǒng)-設(shè)計(jì)材料及要求\data2p.txt','w');</p><p>　　fprintf(fid2,'%12.8f\r\