哈工大信号检测与处理课程报告

1、考核科目:学生所在院（系）：学生所在学科:2017年秋季季学期研穷生课程考核（读书报告、研究报告）信号检测与处理航天学院控制科学与工程学生姓名:学号:学生类别:考核结果17B904012学术型第一部分、信号检测1.相关函数的基础原理相关函数定义为两样本积的数学期望，表示随机信号关联程度、变化程度的量度。是任意样本相应的时间平均值，表示两个样本在不同时间上的相关性。相关函数是信号检测理论中的基础，只有弄清相关性的意义，才能了解后面以相关为基础的一系列方法与原理。特别地，自相关函数定义如下(各态历经下表达式可以由概率平均简化为时间平均如最右表达式):Rxx_1T(T)=Rx(t,t+T)=Ex(t

2、)x(t+T»=1XiX2p(Xi,t;X2,t+T)=¥%,(x(t)x(t+Tyt公式中的期望是在实际中相当于针对时间取的均值，因此相关函数的定义也看作一种对1本身共轲的卷积运算后的平均值：Rxx(T)=；x(tyx(Jo因此，首先讨论卷积的操作与物理意义。卷积物理意义是将信号分解成冲激信号之和，借助系统冲激响应求出系统Nzs对任意激励信号的零状态响应。卷积定义推导如下：将输入信号分解为多个时刻冲激信号的叠加，分别输入并作用于系统如图1。Et,)Nzsjhtt-k)|"Nzsjht-kL.fkit-k、Nzsjfk1：='ht-k；：.n_1n_1二f

3、k-)：、Li：t-kL"：ftNzsJ.fk,)ht-kt"：rtkz0k王则根据以上线性系统输入输出间对应关系可做出如下推导:ft:fati=AfkLt；=；'tkL?；”：t：k1.'：'J£f(小秒一小)一时一5)鼠k丑nX:vfk-1''、：:t-k：.kSn1n-1f(t/工f(k&"6(t-kAT卜r(t怜£f(kSTh(tkS)(2)kWk=e取极限nT,AfTd丁可彳导f(t尸f(T»(t-七XE=f(t)*6(t),即冲激信号与任意输入信号的卷积等于自身；rt(t)=

4、(f)h(t一7)d7=f(t)wh(t),即系统零状态输出等于系统2。冲激响应与输入信号的卷积。也就是卷积运算相当于按照h(t)的分布模式将每个采样点生成的f(t)叠加一遍。卷积应用于LTI(线性时不变系统)，卷积物理意义的图形解释如下图rpCwmitfTltorn-jf弧”啊以门而图2.卷积物理意义的图像解释从以上信号处理中的卷积原理可以看出，关于自身重复事件的信息，体现信号的“记忆能力冲星的频率，或用来估计乐音的音高。自相关也就是信号与自身共轲的卷积可以提供,例如音乐节拍(例如确定节奏)或脉再者从随机变量角度来看相关性，相关函数的引入是为了描述两个随机变量间不同时刻的联系，弥补均值与方差

5、等各阶统计量只描述随机过程的静态数字特征的不足。变量完全线性相关，则会有接近的概率密度分布。以连续随机变量为例，假如随机变量若两个随机x,y且y=5x,那么x,y完全线性相关，X=5的概率和Y=25的概率是相等的，即x,y有相同的概率分布、成线性关系的期望，因此概率角度线性相关性反应的是两个随机过程之间的线性相关程度；从线性空间的角度看相关运算本质是内积运算,而两个向量的内积在线性空间中表示一个向量向另一个向量的投影，因此可以表示两个向量的线性依赖或相似程度,体现出了这种相似程度。求解相关函数Matlab命令如下：dt=.1;t=0:dt:20;x=f(t);y=g(t);a,b=xcorr(

6、x,y,'unbiased');plot(b*dt,a)相关运算以下图3为三种典型自相关函数：a.三次函数;iox106b.正弦函数；c.(0,1)区间均匀分布白噪声.86420-15-10-50510150.20.10-0.1-0.2-0.3-0.4-0.5-20图3.三种典型自相关函数由这些自相关函数可以直观得到以下几点结论：自相关函数通常为偶函数，且Rxx（0尸E<（t力=W十-零点值最大，为信号平均功率;（2）周期序列的自相关函数仍有明显周期性；（3）噪声的自相关性极差，几乎只在原点时刻存在自相关性，与确定信号不同，这也是维纳滤波等相关滤波方

7、法去除噪声的主要原理之一。（左）与瞬态声发射信号卜面图4为将自相关分析运用在平时所处理的钢轨运行过程噪声（右）的结果，可以看出作为一种瞬态信号，声发射信号的相关性同时存在衰减较快与存在一单纯做相关不定周期特性的特点，而轮轨运行产生的噪声也不是传统意义上的白噪声信号。能区分二者，下面选择一种相关性方法与EMD配合使用的方法。-4图4.钢轨运行噪声与声发射信号自相关函数2.基于分段相关EMD的声发射钢轨伤损信号瞬时检测高速情况下分段相关经验模态分解（EMD）的瞬时检测方法重构高速情况下伤损信号的方法，确定伤损信号发生的初始时刻，同时也为提取高速下声发射源特性的研究做出了准备工作。EMD基本理论如下

8、：机械系统的振动信号包含着大量的设备运行的状态信息，对振动信号进行经验模态分解得到的IMF与机械系统的振动模式有着很强的对应关系，所以钢轨振动信号的IMF可以用来分析钢轨伤损。对于每个IMF都应该要满足下面的两个条件：a）在整个数据段内，数据的极值点和数据的过零点的数目必须要相等或者也最多差一个；b）在数据段内的任一点，由数据的局部极大值确定的包络和由数据局部极小值确定的包络的平均值应该为零。为了进一步得到信号的各阶IMF,需要对信号进行经验模态分解，经验模态分解是建立在下面这些假设上：a）信号应该至少有一个极大值和一个极小值；b)通过两个极值点的时间长度来定义特征时间尺度；c)如果信号没有极

9、值点，但是存在一些变形点，可以通过对数据进行求导的方法来获得它的极值点，然后通过积分得到结果。经验模态分解的过程可以分为下面的这些步骤：a)寻找出待分解X(t)的所有极大值数据点Xmax(t)和极小值数据点Xmin(t)；b)分别对所有的极大值数据点和极小值数据点用三次样条曲线进行插值拟合；c)计算出这两条样条曲线的平均值大小m(t)；d)从信号x(t)中减去m(t)h(t)=s(t)-m(t)e)如果h(t)不满足IMF的定义，则还需要进一步分解，把h(t)定义为新的数据x(t),重复上面的步骤，直到h(t)满足IMF的定义，并记c1(t)=h(t),剩下的部分可以表示为r1(t)=x(t)

10、-q(t)f)将ri(t)作为x(t)重复上面的步骤a)到e),直到(t)为单调函数或者足够小的时候，分解过程结束，原始信号可以表示成nx(t)="crni1为了确定这个停止条件，NordenE.Huang等人提出了收敛准则，可以表示为二h1(k二)-"(t)SD=Z2L<名tmh1(kJ)(t)式中：T信号序列的时间长度；;收敛系数，一般可以取0.2到0.3之间的数。传统的小波和滤波等方法是基于固定小波基等，具有局限性，而相关EMD消噪方法以相关滤波及EMD方法为基础，能够适应速度的变化以及相应的条件变化，具有较好的适应性。以静态伤损信号为参考信号，对检测信号分段后

11、进行EMD分解，取分解得到的IMFs分别与参考信号采用互相关函数进行相关性分析，则在声发射有用信号发生的时刻，互相关函数会存在一个极值点，去除不相关的IMF保留相关性强的IMF,最后利用保留下的IMFs重构已经消除噪声的伤损信号，并利用小波时频分析对原始信号与重构信号做出对比。如图12(a)-(b)为18km/h下原始信号及消噪后的重构信号，如图13(a)-(b)为48km/h下原始信号及消噪后的重构信号的小波时频图，从中可以看出去噪前信号噪声及反射信号产生的干扰都很严重，而去噪后可以较完整地保留声发射源产生时刻的不包含反射波形的初始声发射信号。需要说明该方法在噪声信号的成分极多的情况下会出现

12、误检，如图14(a)-(b)为124km/h下原始信号及消噪后的重构信号的小波时频图，其中第一个检测点即为误检。结果表明：说明该方法可以达到去噪并重构声发射信号及确定伤损发生时刻的目的，但其中在噪声极大地情况下存在误检情况。大于阈值的IMF保留，其余舍弃，确定阈值方法：Nref12ref九=1zmax|Rj(E)|R)|,R=(CS广(t)%(t+)dt2NrefT.9Nef附：所有程序见附录图12)ZHKVVeneuaarF)ZHKVVeneuaarFooooO54321oOO521005055215z200图1320050804015030100100402020505010Time(ms

13、)Time(ms)60150143.相关滤波等技术应用(1)相关滤波实例：以下图5为应用相关滤波对包含随机白噪声的正弦波滤波效果:01002003004005006007008009001000输入信号5-5JIIIIIJIII01002003004005006007008009001000实际信号50-5IIIIIIEII01002003004005006007008009001000滤波后信号图5.相关滤波对含随机噪声的正弦波滤波效果相关滤波是借助于参考信号和输入信号中实际信号的相关性，以及参考信号与输入信号中噪声的不相关性，通过求输入信号和参考信号的相关函数，把噪声从输入信号中消除，从而

14、达到滤波效果。原理是利用相关函数在平稳、各态历经下统计平均可用任一样本函数时间平均近似的原理，滤除同频信号，即同频检测原理，此类问题中由于频率已知，最重要的未n知量即相角与幅度%,中0。设原始输出信号为X(t)=X0+£ssin®jt+*)+n(t)包含直流分量、i_0Acos孰t,Asin00t相乘后积分得基波分量、谐波及异波分量、噪声分量。将原始信号分别与Ex,y。直流分量:Exd1T=Eyd、,0X0Asin.0tClt二1limTT二0X0ACOs0tdt：-0基波分量:=0.5SoAcosiR,Eyo=O.5SoAsini'：o噪声分量:Exn_1T=Ey

15、n=ym!n(tJAsin(ctjdt=0直流分量整周期积分后为零，因此随T增大均值周期量会远小于T,故长时间T下可看作积分平均为0;谐波分量与基波分量正交，故积分为0;噪声均值为0,积分为0。而基波积分在时间上平均不受时间影响,为常值。从而由以上关系可用二元方程求解出tan轧=Ey/Ex,5=252。(2)相关滤波实例256图6.相关滤波对含随机噪声及高次谐波的余弦波滤波效果A处的两侧建立测井(3)漏水检测：基于相关算法检测地下管路的漏水点：在怀疑为漏水点XY,分别用仪器测量漏水处漏水的声音，两个测井所测得的漏水声音两条波形相似但交错一段时间。在本实验中，假设XY的距离已知，水速已知，从而水

16、流从X到Y的时间已知。实验中我们用XY两处的相错拍数代替相错时间，取40由于水流的不均匀性，XY两处记录的x(t),y(t)并不恒定，但两者经过40拍后具有很强的相关性。由于二者的距离差在接收信号上会转化为时间差，由相关法测两端信号的互相关可求Ry.=E俨tyt.=E|xtxt-k.当T=k,可知Ry(k)=Ejx2(t=Rxx(0)达到最大，也就是最大值对应拍数为时差，漏水点距离较近观测点S吟。2.8漏水检测图7.漏水检测上图显示了在39步处二者互相关函数最大，与设定的40比较接近。(4)功率谱估计：功率谱描述了信号功率随频率的变化，反

17、映了信号各频率成分的功率分布。2PowerSpectralDensityEstimategoooood224o-euad*7rewBH3口ysnenratceDsrewopFrequency(Hz)yulearPSDestimate50fequency(Hz)DDLyrsnenr,af-ceDsrewop,IZrBorvoneuaetrrewop0000diWelchPowerSpectralDensityEstimate功率谱包括输出自功率谱：Syy(S)=H(W)8x(6)；互功率谱Sxy(6)=H侬)Sxx侬)；第二部分、信号处理4.基于Kalman滤波去除声发射信号采集中的均匀白噪声基

18、于以上分析，由于伤损声发射信号具有瞬发特性，因此信号的相关性不好，用传统的相关滤波效果并不能很好的解决。这里选择在钢材拉伸实验中产生的声发射信号作为原始信号，首先将信号归一化，最大幅值设为1,信号取值范围为卜1,1,并在原始信号的基础上，添加随机产生标准差0.2、方差0.04的均匀分布白噪声，与原始信号叠加并重新归一化，对应操作的简单程序如下：signal1=signal/max(abs(signal);L=2048;a=1;w=(rand(a,L)-0.5)*sqrt(12*0.04);nos_signal=w+signal1;nos_signal1=nos_signal./max(abs(

19、nos_signal);原始信号与添加随机噪声后信号的对比如下图：红色为添加噪声后声发射信号曲线，蓝色为原始信号曲线。图9由于声发射信号具有频带宽（约1MHz）、瞬发性、相关性差等特点，因此下面基于Kalman滤波方法对噪声进行去除。程序如代码5中所示，这里先把设置的噪声方差设为5代入了模型，正常应该单独取出一段噪声信号进行估计。并且去噪后信号幅值通常要小于原始信号，因此对去噪得到信号进行尺度放缩后再与原始信号进行比较如下：红色为去噪后信号，蓝色为原始信号，可以看出滤波的效果超出了原始信号范围，滤波后信号强度弱与原始信号，且有轻微的相位偏移现象，在5MHz采样频率下约6个采样点。图10可求出滤

20、波后信噪比，此时为37.1153,将观测噪声方差改为预设的0.04后滤波效果如下图。图11SNR=10lgPsPn)=20lgVsVNN其中，Vs=maxSignal=0.8661,Vn=RMS(NoisezVn2?=0.035,最后算得信噪比为N64.173,从图中可看出滤波得到的信号与原始信号基本重合，在预估出测量噪声方差后，卡尔曼滤波程序较好的实现了滤波功能。参考文献1知乎.卷积的物理意义是什么EB/OL.http:https:2王翔.功率谱和相关函数在脑血氧监测中的应用.J武汉大学学报.Vol.51.No.3.pp:303307.2005.063宋康.图像自相关函数作为清晰度评价函数的

21、理论分析.J西南交通大学学报.Vol.43.No.5.pp:578582.2008.104 QiushiHao,YiShen,YanWang,XinZhang,DetectionmethodbasedonKalmanfilterforhighspeedraildefectAEsignalonwheel-railrollingrig,NondestructiveTestingandEvaluation.2016,1-17.5 XinZhang,YanWang,KangweiWang,YiShen,HengshanHu.Railcrackdetectionbasedontheadaptivenoi

22、secancellationmethodofEMDathighspeedC.InProc.of2017InstrumentationandMeasurementTechnologyConference(I2MTC),1-6,Turin,Italy,2017.I代码2：%correlation_filtern=0:1024;%产生序列f0=50;%实际信号频率fs=50*1024;%采样频率dt=1/fs;%采样周期T=1/f0;%实际信号周期w=2*pi*f0;%实际信号角频率u=5+5*cos(w*n*dt+30)+6*randn(size(n);%输入信号figure;subplot(3,

23、1,1);plot(u);xlabel('输入信号');realsignal=5*cos(w*n*dt+30);%实际信号subplot(3,1,2);plot(realsignal);xlabel('实际信号');ylabel('幅值');u1=u.*sin(w*n*dt);%滤波后信号v1=u.*cos(w*n*dt);U=sum(u1)*dt*2/T;V=sum(v1)*dt*2/T;A=sqrt(U*U+V*V);p=atan(V/U);%atan(X)isthearctangentoftheelementsofXy=A*cos(w*n*

24、dt+30+p);subplot(3,1,3)plot(y,'r');xlabel('滤波后信号');代码3:%相关滤波例子n=0:511;f0=20;%Hz基频频率N=1;fs=20*512;%Hz采样频率T=1/f0;dt=1/fs;%采样周期w=2*pi*f0;%角频率y=5*cos(w*n*dt)+5*randn(size(n)+5*cos(2*w*n*dt)+5*cos(3*w*n*dt)+5*cos(100*pi*dt*n);%构建信号，初始相位a=0%高次谐波50Hz工频信号figure;plot(y);xlabel('原始信号);real

25、signal=5*cos(w*n*dt);figure;plot(realsignal);xlabel('实际信号');u1=y.*sin(w*n*dt);%原始信号与基准信号相乘v1=y.*cos(w*n*dt);U=sum(u1)*dt*2/(N*T);V=sum(v1)*dt*2/(N*T);amp=sqrt(U*U+V*V);pha=atan(U/V);yout=amp*cos(w*n*dt+pha);%重构得到波形figure;plot(yout);ylabel('幅值');xlabel('滤波后的波形,);代码4:%Water_Leak_De

26、tectionclear;clc;N=1000;n=0:N-1;fs=100;t=n/fs;Lag=100;xx=2+randn(1,length(t);x=xx(100:600);path_X=40;%漏水处距离X为40拍时L=xx(path_X+100:800)-0.5;figure;gridy=xx(1:100),L(100:500);r,lags=xcorr(x,y,Lag,'unbiased');plot(lags,r);xlabel('n');ylabel('rxy(n)');title('漏水检测)代码5:clc;clear

27、;fs=200;n=0:1/fs:1;xn=0.5*cos(2*pi*20*n)+2*cos(2*pi*25*n)+cos(2*pi*50*n)+3*cos(2*pi*80*n)+0.1*randn(size(n);subplot(2,2,1);Hs=spectrum.periodogram;psd(Hs,xn,'fs',fs);subplot(2,2,2);orderl=50;nfft=512;Pxx,f=pcov(xn,orderl,nfft,fs);plot(f,10*log10(Pxx);xlabel('fequency(Hz)');ylabel(

28、9;powerspectraldensity(Db/Hz)');title('pcovPSDestimate');range='half;magunits='db'Pxy,f=pyulear(xn,orderl,nfft,fs,range);subplot(2,2,3);plot(f,10*log10(Pxy);xlabel('fequency(Hz)');ylabel('powerspectraldensity(Db/Hz)');title('yulearPSDestimate');subplot(2,2,4);pburg(xn,orderl,nfft,fs,range);figure;Hs=spectrum.welch;psd(Hs,xn,'