醫(yī)學信號處理實驗報告——兩路信號關(guān)系衡量

1、.電子科技大學生命科學與技術(shù)學院標 準 實 驗 報 告（實驗）課程名稱 生物醫(yī)學信號處理 2018-2019-第2學期電子科技大學教務處制表學生姓名學 號指導教師實驗時間2019.4.16一、實驗室名稱： 品學樓B302二、實驗名稱：兩路信號間的關(guān)系衡量三、實驗學時：2四、實驗原理：1.信噪比(signal-to-noise ratio):是描述信號中有效成分與噪聲成分的比例關(guān)系參數(shù)，單位為dB。假設不含噪聲的信號為sn，sn外加噪聲wn以后的信號為xn=sn+awn，則信號的信噪比定義為snr=10log10(var(xn)var(yn)其中var(xn)代表信號xn的方差在給定信噪比snr

2、的情形下，要求解系數(shù)a，則其計算公式為a=var(xn)varwn10snr102.皮爾遜相關(guān)系數(shù)在統(tǒng)計學中，皮爾遜相關(guān)系數(shù)(Pearson correlation coefficient)，通常用R或表示，是用來度量兩個變量X和Y之間的相互關(guān)系(線性相關(guān))的，取值范圍在 -1,+1 之間。它在學術(shù)研究中被廣泛應用來度量兩個變量線性相關(guān)性的強弱。在作為衡量線性回歸效果時，常使用R2對于隨機變量X和Y，皮爾森相關(guān)系數(shù)的求解公式為： =corX,Y=cov(X,Y)varXvar(Y)其中cov(X,Y) 代表X與Y的協(xié)方差，varX 和var(Y) 分別代表X和Y的方差。當相關(guān)性為 1 時，X與

3、Y的關(guān)系可以表示為Y=aX+b，其中a0；當相關(guān)性為 -1 時，X與Y的關(guān)系可以表示為Y=aX+b，其中a0.1); locs = locs(lo_cs); T_c = mean(diff(locs); plot( xcorr_x(:,i); endend figure(3) plot(0:30,corr_sx,o);disp(實際周期為：);disp(T);disp(相關(guān)測出周期為：);disp(T_c);結(jié)果：實際周期為：200 相關(guān)測出周期為：199.8750圖4 雙頻信號與含有噪聲的雙頻信號圖5 含有/不含有噪聲的信號線性相關(guān)結(jié)果（有偏）圖6 相關(guān)系數(shù)與系數(shù)a關(guān)系圖（三）程序：clea

4、r all;N=1000; %長度Fs=50; %采樣頻率n=0:N-1; t=n/Fs; %時間序列A=0.4;A1=0.5;A2=0.6; %衰減系數(shù)c0=340; %c0d1=620d2=500 %請自己給兩個距離的參數(shù)，不要太小即可t1=d1/c0;t2=(d1+2*d2)/c0;tc=2*(d1+d2)/c0;Lag=500; %最大延遲樣點數(shù)pt=sinc(2*pi*t); %原信號xt = pt + A1*sinc(2*pi*(t-tc) ;figure(1)plot(1:N)/Fs,xt);yt = A*sinc(2*pi*(t-t1)+ A2*sinc(2*pi*(t-t2)

5、;figure(2)plot(1:N)/Fs,yt);Rpp,lags = xcorr(pt,biased); %p(t)自相關(guān)Rxx,lagx = xcorr(xt,biased); %x(t)自相關(guān)Ryy,lagt = xcorr(yt,biased); %y(t)自相關(guān)Rxy,lagy = xcorr(xt,yt,biased); %x(t)與y(t)互相關(guān)rts=lags/Fs;rtx=lagx/Fs;rtt=lagt/Fs;rty=lagy/Fs;figure(3)pks,locs=findpeaks(Rxy);lo_cs = find(pks1E-3);locs = (locs(l

6、o_cs)-999)/Fs;pks = pks(lo_cs);plot(locs,pks,o);hold on;plot(rty,Rxy);figure(4)subplot(2,2,1)plot(rts,Rpp)title(%p(t)自相關(guān))xlabel(時間/s)subplot(2,2,2)plot(rtx,Rxx)title(x(t)自相關(guān))xlabel(時間/s)subplot(2,2,3)plot(rtt,Ryy)title(%y(t)自相關(guān))xlabel(時間/s)subplot(2,2,4)plot(rty,Rxy)title(%x(t)與y(t)互相關(guān))xlabel(時間/s)d

7、1_c = (locs(3) * c0 - locs(2) * c0)/2d2_c = (locs(4) * c0 - d1_c)/2 + (-locs(1) * c0 - d1_c)/2)/2結(jié)果：圖7 三個信號自相關(guān)及x，y互相關(guān)圖d1d2d1（測）d2（測）620500612499.8650420646418.2800700792.2702.1420650411.4651.1表1 d1，d2測量情況表（四）程序：clear all ;close=load(eegclose);open=load(eegopen);close1 = close.eegclose(:,8*2);close2

8、= close.eegclose(:,9*2);open1 = open.eegopen(:,8*2);open2 = open.eegopen(:,9*2);fs = 1000;t = -5:1/fs:5;s1 = sin(50*2*pi*t) + 2*sin(250*2*pi*t) + randn(1,length(t);s2 = sin(150*2*pi*t) + 4*sin(250*2*pi*t) + randn(1,length(t);figure(1)pwelch(s1,hamming(512),50,1024,fs);figure(2)pwelch(s2,hamming(512)

9、,50,1024,fs);figure(3)mscohere(s1,s2,hamming(512),50,1024,fs);% 在功率譜中標注出(1-3Hz)、(4-7Hz)、(8-13Hz)、(14-30Hz)locate_mid = 1.5,5.5,10.5,22,50,100;fs2 = 250;n = 128;nfft = 512;noverlap = 1/2*n;figure(4)Cxy,F = mscohere(open1,open2,blackman(n),noverlap,nfft,fs2);plot(F,Cxy);title(不同路睜眼信號頻域相干)set(gca,XTick

10、mode,manual,Xtick,locate_mid)figure(5)Cxy,F = mscohere(close1,close2,blackman(n),noverlap,nfft,fs2);plot(F,Cxy);title(不同路睜眼信號頻域相干)set(gca,XTickmode,manual,Xtick,locate_mid)figure(6)Cxy,F = mscohere(close1,open1,blackman(512),256,1024,fs2);plot(F,Cxy);title(同路信號頻域相干)set(gca,XTickmode,manual,Xtick,loc

11、ate_mid)結(jié)果：圖8 信號1的功率譜圖9 信號2的功率譜圖10 兩信號的頻率相干圖圖11不同路睜眼信號的頻率相干圖圖12不同路閉眼信號的頻率相干圖 圖13同路信號的頻率相干圖(1-3Hz)(4-7Hz)(8-13Hz)(14-30Hz)35Hz-37Hz均值0.0744381080.0302513160.0301312270.0683740770.087902008最大值0.3775762960.0884963430.0929708640.3114252990.319218072表2 同路信號相干結(jié)果表（五）程序：clear all;clca = 1 2 5 6 7;b = 2 5 6;

12、s1 = Icorr(a,b)s2 = xcorr(a,b)function s = Icorr(s1,s2)% 計算兩個函數(shù)的線性相關(guān)s1_s = size(s1);s2_s = size(s2);if s1_s(1) = 1 s1 = s1;endif s2_s(1) = 1 s2 = s2;ends1_len = length(s1);s2_len = length(s2);s_len = max(s1_len,s2_len);dlen = abs(s1_len-s2_len);if s_len = s1_len; s2 = s2,zeros(1,dlen);else s1 = s1,z

13、eros(1,dlen);endy = zeros(s_len,s_len*2-1);for i = 1:s_len y(i,s_len+i-1:-1:i) = s1;ends = s2*y;s = s(1:s_len*2-1-dlen);結(jié)果：Icorr： 14 47 82 65 42 17 6Xcorr：0.0000 0.0000 6.0000 17.0000 42.0000 65.0000 82.0000 47.0000 14.000010、 總結(jié)及心得體會：（一）1. weight值隨著信噪比減小而升高。2. 構(gòu)造信號（x=s+weight*n），信噪比越高相關(guān)系數(shù)越大。在20dB時已

14、經(jīng)接近最大為1，在-20dB時接近最小值0.04。構(gòu)造信號（x=weight*n-s）時，在-20dB時已經(jīng)接近最小為-1，在20dB時接近最大值-0.04。這說明，如果要從含噪聲的信號中檢測到信號，一般手段的極限在-20,20dB之間。3.相關(guān)系數(shù)的符號代表了進行相關(guān)計算兩信號序列的變化趨勢關(guān)系。如為正，說明一個信號序列中增大時，另一個相應的增大，反之減小，稱為具有正相關(guān)性。符號為負則相反，稱為具有負相關(guān)性。 （二）1.將信號與含噪聲的信號進行線性相關(guān)可以有效地識別出周期。設想該性質(zhì)是由于加入的是randn函數(shù)生成的白噪聲信號。改為加入rand函數(shù)生成的噪聲，自相關(guān)的方法仍能準確體現(xiàn)出信號的

15、周期性。變化的只有相關(guān)系數(shù)的變化趨勢變得趨于線性，不過這一點在意料之中。2.如果使用無偏的自相關(guān)，會收獲更加明顯的周期性和整齊的峰值，但結(jié)果差別不大。（三）用延時相關(guān)估計距離是可行的，誤差并不大。d2的誤差規(guī)律不明顯，但是d1的測量值總是較d1的真實值要小，該方法有可能存在系統(tǒng)誤差。（四）1.兩信號頻率相干圖的峰值出現(xiàn)的位置也是兩信號功率譜均具有峰值的頻率。2. 為了方便分析不同路信號相干圖像，程序中選取128點為一段，并使用了主瓣更寬的blackman使圖像更加平滑。1）不同路同類型信號的相干圖中可以觀察到在功率譜中標注出(1-3Hz)、(4-7Hz)、(8-13Hz)、(14-30Hz)的范圍均有對應峰值。36Hz左右有未知原因的峰值， 46Hz后應該都是干擾性的噪聲所以相干值保持很高。這說明這些位置對應的頻率分量大量包含在信號中。2）可以看出與睜眼信號相比，閉眼信號除了波(1-3Hz)的值較低以外，