機(jī)械工程測試技術(shù)基礎(chǔ)

1、學(xué)習(xí)要求:,1 數(shù)字信號處理概述 2 了解信號模數(shù)轉(zhuǎn)換原理　　 3 掌握信號采樣定理，能正確選擇采樣頻率 4 了解信號截斷、能量泄露、柵欄效應(yīng)等現(xiàn)象　 5 DFT與FFT 6 柵欄效應(yīng)與窗函數(shù)7 常用的數(shù)字信號處理算法,1 數(shù)字信號處理概述,1)數(shù)字信號處理的主要研究內(nèi)容,數(shù)字信號處理主要研究用數(shù)字序列來表示測試信號，并用數(shù)學(xué)公式和運(yùn)算來對這些數(shù)字序列進(jìn)行處理。其主要內(nèi)容包括數(shù)字波形分析、幅值分析、頻譜分析和數(shù)字濾波。,2)

2、測試信號數(shù)字化處理的基本步驟,3) 數(shù)字信號處理的優(yōu)勢,(1)用數(shù)學(xué)計算和計算機(jī)顯示代替復(fù)雜的電路和機(jī)械結(jié)構(gòu),(2)計算機(jī)軟硬件技術(shù)發(fā)展的有力推動,a)多種多樣的工業(yè)用計算機(jī)。,b)靈活、方便的計算機(jī)虛擬儀器開發(fā)系統(tǒng),2 模數(shù)(A/D),1) A/D轉(zhuǎn)換過程,采樣――利用采樣脈沖序列p(t)，從連續(xù)時間信號x(t)中抽取一系列離散樣值，使之成為采樣信號x(nTs)的過程．,編碼――將離散幅值經(jīng)過量化以后變?yōu)槎M(jìn)制數(shù)字的過程。,如4位碼，

3、只表示24=16種不同的信號幅度，這些幅度稱為量化電平。當(dāng)離散時間信號幅度與量化電平不相同時，就要以最接近的一個量化電平來近似它。所以經(jīng)過A/D變換器后，不但時間離散化了，而且幅度也量化了，產(chǎn)生一個二進(jìn)制流。,2) A/D轉(zhuǎn)換器的技術(shù)指標(biāo),分辨率 用輸出二進(jìn)制數(shù)碼的位數(shù)表示。位數(shù)越多，量化誤差越 小，分辨力越高。常用有8位、10位、12位、16位等。 轉(zhuǎn)換精度 轉(zhuǎn)換精度?LSB/2 LSB: Le

4、ast Significant Bit 最低有效位轉(zhuǎn)換速度 完成1次轉(zhuǎn)換所用的時間,如100ms(10Hz);10us(100kHz) 模擬信號的輸入范圍;如5V， +/-5V，10V，+/-10V等。,3 采樣定理,1) 采樣信號的頻譜,采樣過程是將采樣脈沖序列p(t)與信號x(t)相乘來．,2) 頻混現(xiàn)象,頻域解釋,,3)采樣定理,為保證采樣后信號能真實(shí)地保留原始模擬信號信息，信號采樣頻率必須至少為原信號中最高頻率成分的2倍。

5、這是采樣的基本法則，稱為采樣定理。,fs＞2fmax,需要注意的是，在對信號進(jìn)行采樣時，滿足了采樣定理，只能保證不發(fā)生頻率混疊，只能保證對信號的頻譜作逆傅立葉變換時，可以完全變換為原時域采樣信號xs(t)，而不能保證此時的采樣信號能真實(shí)地反映原信號x(t)。工程實(shí)際中采樣頻率通常大于信號中最高頻率成分的3到5倍。,頻混現(xiàn)象實(shí)驗(yàn)：,頻混計算：,正常,混迭頻率=Fs-信號頻率,A/D采樣前的抗混迭濾波：,計算機(jī),,顯示,4 信號的截斷、能

6、量泄漏,為便于數(shù)學(xué)處理，通常對截斷的信號做周期延拓，得到虛擬的無限長的信號。,當(dāng)運(yùn)用計算機(jī)進(jìn)行測試信號處理時，不可能對無限長的信號進(jìn)行測量和運(yùn)算，而是取其有限的時間片段進(jìn)行分析，這個信號截取過程成為信號的截斷。,周期延拓后的信號與真實(shí)信號是不同的，下面我們就從數(shù)學(xué)的角度來看這種處理帶來的誤差情況。,設(shè)有余弦信號x(t), 用矩形窗函數(shù)w(t)與其相乘，得到截斷信號: y(t) =x(t)w(t),將截斷信號譜 XT(ω)與原始信號譜

7、X(ω)相比較可知，它已不是原來的兩條譜線，而是兩段振蕩的連續(xù)譜. 原來集中在f0處的能量被分散到兩個較寬的頻帶中去了，這種現(xiàn)象稱之為頻譜能量泄漏。,能量泄漏實(shí)驗(yàn)：,6.5 DFT與FFT,1)離散傅立葉變換,離散傅里葉變換（Discrete Fourier Transform）一詞是為適應(yīng)計算機(jī)作傅里葉變換運(yùn)算而引出的一個專用名詞。,周期信號xT(t)的傅里葉變換：,對周期信號xT(t)采樣，將離散序列xT(n),將積分轉(zhuǎn)為集合：,

8、按上式，用計算機(jī)編程很容易計算出指定頻率點(diǎn)的值：,f=？ //計算的頻率點(diǎn)Fs=5120N=1024dt=1.0/Fspi=3.1415926 XR=0XI=0,For n=0 To N-1 XR=XR+x(n)*cos(2*pi*f*n*dt)*dt XI=XI+x(n)*sin(2*pi*f*n*dt)*dtNextA=sqr(XR*XR+XI*XI) Q=atn(XI/XR),VBScript 樣例

9、,展開，得連續(xù)傅立葉變換計算公式：,連續(xù)傅立葉變換編程計算實(shí)驗(yàn)：,采樣信號頻譜是一個連續(xù)頻譜，不可能計算出所有頻率點(diǎn)值，X( f )只能離散取值，頻率取樣間隔定義為：,Δf=fs/N,頻率取樣點(diǎn)為｛0, Δf, 2Δf, 3Δf, ....｝，有：,2) 快速傅立葉變換 (FFT),FFT是實(shí)施DFT的一種有效算法，通過仔細(xì)選擇和重新排列中間結(jié)果，速度上較DFT有明顯的優(yōu)點(diǎn)。,展開各點(diǎn)的DFT計算公式：,XR(1)=x(0).cos(2

10、pi*0*1/N)+x(1).cos(2pi*1*1/N)+x(2).cos(2pi*2*1/N)…..XR(2)=x(0).cos(2pi*0*2/N)+x(1).cos(2pi*1*2/N)+x(2).cos(2pi*2*2 /N)…..,其中有大量cos、sin項的重復(fù)計算，F(xiàn)FT就是用技巧減少這些重復(fù)計算。,當(dāng)采樣點(diǎn)數(shù)為1024點(diǎn),DFT要求一百萬次以上計算量，而FFT則只要求10240次。,FFT的Matlab實(shí)現(xiàn),fs=1

11、000 t=0:1/fs:0.6; f1=100; f2=300; x=sin(2*pi*f1*t)+sin(2*pi*f2*t); subplot(711) plot(x); title(‘f1(100Hz)\f2(300Hz)的正弦信號，初相0’) xlabel(‘序列（n）’) grid on number=512 y=fft(x,number); k=0:length(y)-1; f=fs*k/leng

12、th(y); subplot(713) plot(f,abs(y)); title('f1\f2的正弦信號的FFT（512點(diǎn)）') xlabel('頻率Hz') grid on,x=x+randn(1,length(x)); subplot(715) plot(x); title('原f1\f2的正弦信號(含隨機(jī)噪聲) ') xlabel(‘序列（n）’) grid

13、ony=fft(x,number); k=0:length(y)-1; f=fs*k/length(y); subplot(717) plot(f,abs(y)); title('原f1\f2的正弦信號(含隨機(jī)噪聲)的FFT(512點(diǎn))') xlabel('頻率Hz') grid on,6 柵欄效應(yīng)與窗函數(shù),1)柵欄效應(yīng),采樣信號的頻譜，常采用FFT算法進(jìn)行計算，設(shè)數(shù)據(jù)點(diǎn)數(shù)為：,N = T

14、/dt = T. fs,則計算得到的離散頻率點(diǎn)為:,Xs(fi) ， fi = k.fs/N , k = 0，1，2，.....，N/2,相當(dāng)于透過柵欄觀賞風(fēng)景，只能看到頻譜的一部分，而其它頻率點(diǎn)看不見，稱此種現(xiàn)象為柵欄效應(yīng)。,如果信號中的頻率分量與頻率取樣點(diǎn)不重合，則只能按四舍五入的原則，取相鄰的頻率取樣點(diǎn)譜線值代替。,柵欄效應(yīng)誤差實(shí)驗(yàn)：,2) 能量泄漏與柵欄效應(yīng)的關(guān)系,頻譜的離散造成了柵欄效應(yīng)，譜峰越尖銳，產(chǎn)生誤差的可能性就

15、越大。,例如，余弦信號的頻譜為線譜。信號頻率與頻譜離散取樣點(diǎn)不等時，柵欄效應(yīng)的誤差為無窮大。,,因信號截斷產(chǎn)生了能量泄漏，即使信號頻率與頻譜離散取樣點(diǎn)不等，也能得到一個該頻率分量的近似值，因此能量泄漏誤差不完全有害的。如沒有信號時域截斷產(chǎn)生的能量泄漏誤差，頻譜離散取樣造成的柵欄效應(yīng)誤差不能接受。,能量泄漏分主瓣泄漏和旁瓣泄漏，主瓣泄漏可以減小因柵欄效應(yīng)帶來的譜峰幅值估計誤差，有其好的一面，而旁瓣泄漏則是完全有害的。,3) 常用的窗函數(shù),

16、(1) 矩形窗,(2) 三角窗,(3) 漢寧窗,窗函數(shù)在減小柵欄效應(yīng)誤差中的作用實(shí)驗(yàn)：,總結(jié)：,信號截斷?能量泄漏,FFT?柵欄效應(yīng),從克服柵欄效應(yīng)誤差的角度看，能量泄漏是有利的。,通過加窗加大能量泄漏，減小柵欄效應(yīng)誤差：,7 常用的數(shù)字信號處理算法,數(shù)字信號處理是測試技術(shù)中最常用和最需要應(yīng)用工程師掌握的部分，無論開發(fā)簡單或復(fù)雜的測控系統(tǒng)或分析儀器，都會用到數(shù)字信號處理知識。,,,,,信號時域波形分析是最常用的信號分析手段，用示波器、萬

17、用表等普通儀器顯示信號波形就可以特征參數(shù)。,1)時域波形參數(shù)計算,,波形分析的應(yīng)用,超門限報警,① 峰值P，雙峰值Pp-p,② 均值,U=0For K = 0 To N U=U+data(k)NextU=U/N,③ 均方值,E2=0For K = 0 To NE2=E2+data(k)*data(k)NextRMS=sqr(E2/N),④方差,V2=0For K = 0 To N V2=V2+(data(k)

18、-U)*(data(k)-U)NextV=V2/N,⑤周期T,n=0AT=0.8*PFor K = 2 To N If data(k-1)AT And data(k+2)>AT Then ti(n)=K n=n+1 End IfNextT=(ti(2)-ti(1))*dt,兩個序列的卷積：c(n)=a(n)*b(n)= 如a和b的腳標(biāo)范圍為：[a1,a2]與[

19、b1,b2]，則c的腳標(biāo)范圍為[a1+b1,a2+ b2] ，但matlab的索引是從1開始的，設(shè)a和b都是從n=1開始, 知道最后的范圍即可.function c=convolution(a,b)M=length(a); %a序列的長度N=length(b); %b序列的長度 c=zeros(1,N+M-1); %用來存放卷積的結(jié)果for n=2:(M+N) %卷積之后的索引本應(yīng)該的范圍 tmp_

20、max=min(n-1,M); %根據(jù)求和表達(dá)式，由1=<k<=M 且 1=<n-k<=N tmp_min=max(n-N,1); %確定最終求和變量k的取值 for k =tmp_min:1:tmp_max; c(n-1) = c(n-1) +a(k)*b(n-k); %求和，同時由于c從1開始索引， end

21、 %將n的范圍2：M+N對應(yīng)到c的1：M+N-1中去end,(2) 卷積的計算,(2) 卷積的計算,調(diào)用matlab自帶卷積函數(shù)conv驗(yàn)證：>> x = 1:9;>> y=2:8;>> conv(x,y)ans = 2 7 16 30 50 77 112 147 182 197 200

22、 190 166 127 72>> convolution(x,y)ans = 2 7 16 30 50 77 112 147 182 197 200 190 166 127 72,(3) 數(shù)字相關(guān)函數(shù)計算,變量之間的相依關(guān)系稱為相關(guān)。信號之間的相似關(guān)系稱為相關(guān)函數(shù)。,相關(guān)分析的工程應(yīng)用,案例：機(jī)械加工表面粗糙度自相關(guān)分析,

23、被測工件,相關(guān)分析,?提取出回轉(zhuǎn)誤差等周期性的故障源。,案例：地下輸油管道漏損位置的探測,案例： AGV小車定位，聲位筆定位,實(shí)驗(yàn)：自相關(guān)分析,,,如果a和b的腳標(biāo)范圍為：與，則c的腳標(biāo)范圍為，但是由于matlab的索引是從1開始的，這里假設(shè)a和b都是從n=1開始。知道最后的范圍即可：function c =correlation(a,b)M = length(a);N = length(b);c=zeros(1,M+N-1);

24、 %用來存放卷積的結(jié)果for n=(-N+1):M-1 %相關(guān)運(yùn)算本應(yīng)該的范圍 tmp_min = max(1,1-n); %根據(jù)求和表達(dá)式，由1=> xcorr(x,y)=-0.0 -0.0 8.0 23.0 44.0 70.0 100.0 133.0 168.0 203.0 238.0 193.0 150.0 110.0 74.0 43.0

25、18.0>> correlation(x,y) =8 23 44 70 100 133 168 203 238 193 150 110 74 43 18,計算公式：,dt=0.1;t=[0:dt:100];x=cos(t);[a,b]=xcorr(x,‘unbiased’);plot(b*dt,a) 對于互相關(guān)函數(shù)，把[a,b]=

26、xcorr(x,'unbiased');改為 [a,b]=xcorr(x,y,'unbiased');便可。,在matlab中實(shí)現(xiàn)自相關(guān)的無偏估計,延遲序數(shù),(4) 功率譜估計與Matlab,功率譜：隨機(jī)信號的功率譜反映的是隨機(jī)信號的頻率成分及各成分的相對強(qiáng)弱。功率譜估計：基于有限的數(shù)據(jù)尋找信號、隨機(jī)過程或系統(tǒng)的頻率成分。分參數(shù)化和非參數(shù)化2大類方法。非參數(shù)化方法又稱經(jīng)典譜估計，根據(jù)截取的N個

27、樣本數(shù)據(jù)估計出其功率譜．如周期圖法、自相關(guān)法等，主要缺陷是描述功率譜波動的方差性能較差，頻率分辨率低;參數(shù)化譜估計又稱現(xiàn)代譜估計，如AR模型法、移動平均模型法(MA模型法)、自回歸移動平均模型法(ARMA模型法)等．,根據(jù)維納—辛欽定理，對于離散隨機(jī)信號有: Rx(m)為離散隨機(jī)信號自相關(guān)函數(shù)，Sx(ej?)為功率譜密度．如獲得Rx(m) ，Sx(ej?)就不難求出, 故該法的實(shí)質(zhì)是相關(guān)函數(shù)估值的傅立葉變換即為功率譜密度．平穩(wěn)

28、隨機(jī)信號的相關(guān)函數(shù)是確定性函數(shù)，其功率譜也是確定的． 由平穩(wěn)隨機(jī)離散信號的有限個離散值x(0)，x(1)，…，x(N－1)求出自相關(guān)函數(shù)．然后在(-M，M)內(nèi)對Rx(m)作Fourier變換，得到功率譜,用相關(guān)函數(shù)獲得功率譜(BT法),功率譜可用FFT算法求出．M趨于N時，上式求得的Rx(m)方差較大，可靠性差，譜估計質(zhì)量下降．為此Rx(m)上加窗函數(shù)w(m) 使方差大的自相關(guān)函數(shù)估值加權(quán)小，以減少對譜估計的影響．

29、,功率譜是非負(fù)值，加窗后的功率譜也是非負(fù)值．因此窗函數(shù)的傅立葉變換必須是非負(fù)值，這樣窗函數(shù)必須是偶對稱的，上式中的三角形窗是常被采用． 該法的缺點(diǎn)是M趨于N時，R(m)的方差很大，使譜估計質(zhì)量下降；由R(m)得到的S(w)不一定為正值，從而可能失去功率譜的物理意義．該法對由序列x(n)估計出的自相關(guān)函數(shù)R(m)進(jìn)行傅立葉變換，得到x(n)的功率譜估計．當(dāng)延遲與數(shù)據(jù)長度相比很小時，有良好的估計精度,Fs=500; n=0:1/Fs:1

30、; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+0.1*randn(size(n)); window=boxcar(length(xn));nfft=512;cxn=xcorr(xn, 'unbiased');CXk=fft(cxn, nfft);Pxx=abs(CXk);index=0:round(nfft/2－1);k=index*Fs/nfft; plot_Pxx=10*

31、log10(Pxx(index+1));figure(3);plot(k, plot_Pxx);,基于Matlab實(shí)現(xiàn)的程序:,Schuster首先提出的根據(jù)各態(tài)歷經(jīng)隨機(jī)過程功率譜的定義進(jìn)行的譜估計, 把隨機(jī)序列x(n)的N個觀測數(shù)據(jù)視為能量有限的序列，直接計算x(n)的DFT，得x(k)，然后取幅值的平方并除以N，作為序列x(n)真實(shí)功率譜的估計．由平穩(wěn)隨機(jī)信號X(n)的有限個觀察值x(0)，x(1)，…，x(N-1)求其傅立葉變

32、換然后進(jìn)行譜估計 ，M=N－1時，相關(guān)函數(shù)法與周期圖法估計出的功率譜相同；M＜N－1時，相關(guān)函數(shù)法的偏差大于周期圖法,窗函數(shù)滿足一定條件時是漸進(jìn)無偏估計；其方差小于周期圖法；分辨率比周期圖法低，與窗函數(shù)的選擇有關(guān)．可見，窗函數(shù)的選擇對譜估計的修正程度有影響．,周期圖法(Periodogram),周期圖法因與序列的頻譜有直接對應(yīng)關(guān)系，且可采用FFT來計算而被廣泛應(yīng)用．

33、但需對無限長的平穩(wěn)隨機(jī)序列進(jìn)行截斷，相當(dāng)于加矩形窗而成為有限長數(shù)據(jù)，這也意味著對自相關(guān)函數(shù)加三角窗，使功率譜與窗函數(shù)卷積，產(chǎn)生頻譜泄漏，容易使弱信號的主瓣被強(qiáng)信號的旁瓣所淹沒，造成頻譜的模糊和失真，較低譜分辨率．,Nfft=512,Nfft=1024,Fs=600;n=0∶1/Fs∶1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+0．1*randn(size(n)); window=boxcar(len

34、gth(xn));nfft=512; ［Pxx，f］=periodogram(xn，window，nfft，F(xiàn)s);figure(1); plot(f，10*log10(Pxx));window=boxcar(length(xn)); nfft=1024; ［Pxx，f］=periodogram(xn，window，nfft，F(xiàn)s);figure(2); plot(f，10+log10(Pxx));,周期圖法基于Mat

35、lab實(shí)現(xiàn)的程序:,周期圖法不滿足一致估計條件，需改進(jìn)。改進(jìn)目的主要是將周期圖進(jìn)行平滑，減小估計的協(xié)方差，得到一致譜估計。數(shù)據(jù)長度N太大時，周期圖的功率譜估計的譜曲線起伏加?。蝗鬘太小，譜的分辨率又不好。改進(jìn)的主要方法有 Bartlett平均周期圖法 將N點(diǎn)的有限長序列x(n)分段求周期圖再平均．將長度為N的數(shù)據(jù)分為L段，每段長度為M．先對每段數(shù)據(jù)用周期圖法進(jìn)行譜估計，然后對L段求平均得到長度為N的數(shù)據(jù)的功率譜．可得功

36、率譜為 該估計的分辨率 式中由于M遠(yuǎn)小于N，故其分辨率比周期圖法低．可見，該法方差的改善是以犧牲分辨率為代價的．,經(jīng)典譜估計的改進(jìn),Welch平滑平均周期圖法 其改進(jìn)了Bartlett法方差性能和分辨率．基本原理是先對隨機(jī)序列分段，使每一段有部分重疊，然后對每一段數(shù)據(jù)用一個合適的窗函數(shù)進(jìn)行平滑處理，最后對各段譜求平均．得功率譜其中 ，w(n)是窗函數(shù)，其為矩形窗、重疊5

37、0%時，該估計的分辨率 Welch法允許各段數(shù)據(jù)交疊，數(shù)據(jù)段數(shù)L會增加，使方差得到更大的改善，但數(shù)據(jù)的交疊增加了各段數(shù)據(jù)的相關(guān)性，各段數(shù)據(jù)不再是完全獨(dú)立的了,方差的減小不會達(dá)到理論程度．另外采用合適的窗函數(shù)可以減少信號的頻譜泄露，同時也可以增加譜峰的寬度，從而提高分辨率．,Fs=600;n=0∶1/Fs∶1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));

38、nfft=512; window=boxcar(100); %矩形窗windowl=hamming(100); %海明窗window2=blackman(100);%blackman窗 noverlap=20; %數(shù)據(jù)重疊 range=‘half’; %頻率區(qū)間為［0 Fs/2］，計算一半的頻率 ［Pxx，f］=pwelch(xn，

39、window，noverlap，nfft，F(xiàn)s，range); ［Pxxl，f］=pwelch(xn，windowl，noverlap，nfft，F(xiàn)s，range); ［Pxx2，f］=pwelch(xn，window2，noverlap，nfft，F(xiàn)s，range); plot_Pxx=10*log10(Pxx);plot_Pxxl=10*log10(Pxxl);plot_Pxx2=10*log10(Pxx2);figur

40、e(1)plot(f，plot_Pxx);figure(2)plot(f，plot_Pxxl);figure(3) plot(f，plot_Pxx2);,基于MATLAB實(shí)現(xiàn)的程序?yàn)?,矩形窗,漢明窗,布萊克曼窗,(5)數(shù)字濾波,數(shù)字濾波器是利用離散時間系統(tǒng)特性對輸入信號進(jìn)行加工處理，把輸入序列x(n) 變換成一定的輸出序列 y(n)，從而達(dá)到改變信號頻率構(gòu)成的目的。,y(n)=x(n)*h(n),y(k)=h0x(k)+h1

41、x(k+1)+h2x(k+2)+…+ hmx(k+m),k=0,1,........,設(shè)數(shù)字濾波器的脈沖響應(yīng)序列為{h0, h1,h2,…,hm}, 則：,數(shù)字濾波器系數(shù)有專用軟件計算，這里我們提供一個網(wǎng)絡(luò)計算工具：,樣例信號：,,將數(shù)字濾波器系數(shù)帶入公式中就可以對信號進(jìn)行數(shù)字濾波。,For i=0 To N-7 Y(i)=-0.064X(0)+0.041X(1)+0.301X(2)+0.454X(3) +0.30