1、生物醫(yī)學信號處理評分大理大學實驗報告 課程名稱 生物醫(yī)學信號處理 實驗名稱 數(shù)字相關和數(shù)字卷積 專業(yè)班級 姓 名 羽卒蘭cl 學 號 實驗日期 實驗地點 20152016學年度第 3 學期一、 實驗目的熟悉數(shù)字相關的運算，初步在信號處理中應用相關技術。二、實驗環(huán)境 1、硬件配置：Intel（R） Core(TM) i5-4210U CPU 1.7GHz 1.7GHz 安裝內存（RAM）：4.00GB 系統(tǒng)類型：64位操作系統(tǒng) 2、軟件環(huán)境：MATLAB R2013b軟件三、實驗內容(包括本實驗要完成的實驗問題及需要的相關知識簡單概述)已知發(fā)射波形，利用相關技術，在有強背景噪聲的情況下檢測回波的

2、延時和強度。首先使用已知信號模版及其若干次衰減延遲生成仿真回波波形，然后與白噪聲背景疊加，構造仿真信號。然后計算模版與仿真信號的相關函數(shù)，判斷回波位置及相對強度。四、實驗結果與分析 （包括實驗原理、數(shù)據(jù)的準備、運行過程分析、源程序（代碼）、圖形圖象界面等）實驗原理相關可以從時域角度表現(xiàn)信號間的相似（關聯(lián)）程度，是統(tǒng)計信號處理最基本的手段之一。設有離散信號x(n)和y(n)，線性相關函數(shù)定義為 實際采集的信號總是有限長度，用有限的樣本估計相關（自相關）函數(shù) 求和項總數(shù)不是N而是N-|m|，因為當n=N-|m|-1時，n+|m|=N-1。此時xn+m已經(jīng)到了數(shù)據(jù)邊沿。這種估計是漸進無偏估計和一致估

3、計。計算中，只要將其中一個序列反轉，就可以用計算線性卷積的程序計算線性相關 因此可以用FFT來加速相關運算，即對序列補零后，用循環(huán)相關計算線形相關，然后用循環(huán)卷積的快速算法計算循環(huán)相關，得到最終結果。 源代碼：% 估計兩個相似信號間的時間延遲np = 0:99; %從0到99取值% p = sin(pi/5*np); % 正弦% p = exp(-0.06*np); % 指數(shù)衰減% p = sin(pi/5*np).*exp(-0.06*np); % 指數(shù)衰減正弦p = ones(size(np); % 方波figure;subplot(2,2,1); plot(np,p);%數(shù)據(jù)圖n = 1

4、:1000; %定義1000點長的隨機信號w = randn(size(n); %產(chǎn)生1000點長的隨機信號噪聲s = zeros(size(n); %產(chǎn)生1000點長的有用信號%通過改變衰減系數(shù)對已知信號模版進行若干次衰減，生成仿真回波A = 3; % 衰減系數(shù)s(100:199) = s(100:199)+A*p;s(500:599) = s(500:599)+(A/3)*p;s(800:899) = s(800:899)+(A/3/3)*p;x = s+w; %仿真回波與白噪聲疊加構造仿真信號figure;subplot(3,1,1); plot(n,w); title(Noise);

5、%作出噪聲信號的圖，并添加標題subplot(3,1,2); plot(n,s); title(Signal); %作出有用信號的圖，并添加標題subplot(3,1,3); plot(n,x); title(Signal with Noise); %作出觀測信號圖，并添加標題p = p,zeros(1,length(x)-length(p); %如果要求歸一化相關系數(shù)（相干系數(shù)），兩個序列要同樣長Rps = xcorr(s,p,coeff); %對有用信號做互相關函數(shù)Rpw = xcorr(w,p,coeff); %對噪聲信號做互相關函數(shù)Rpx = xcorr(x,p,coeff); %對觀

6、測信號做互相關函數(shù)n2 = (n(1)-np(end):(np(end)-n(1); %線性相關的范圍figure;subplot(3,1,1); plot(Rps); title(Rpw of p(n) and s(n); %作出有用信號線性相關后的圖，并添加標題subplot(3,1,2); plot(Rpw); title(Rps of p(n) and w(n); %作出噪聲信號線性相關后的圖，并添加標題subplot(3,1,3); plot(Rpx); title(Rpx of p(n) and x(n); %作出觀測信號線性相關后的圖，并添加標題結果顯示：1.模板為方波信號，A=

7、3，噪聲均值為0，方差為1a） b）c） 圖1 方波信號的數(shù)據(jù)圖a、信號圖b和線性相關圖c2. 模板為正弦信號，A=3，噪聲均值為0，方差為1a） b）c） 圖2 正弦信號的數(shù)據(jù)圖a、信號圖b和線性相關圖c3. 模板為指數(shù)衰減信號，A=3，噪聲均值為0，方差為1a） b） c） 圖3 指數(shù)衰減信號的數(shù)據(jù)圖a、信號圖b和線性相關圖c4. 模板為指數(shù)衰減正弦信號，A=3，噪聲均值為0，方差為1a） b） c） 圖4 指數(shù)衰減正弦信號的數(shù)據(jù)圖a、信號圖b和線性相關圖c分析： 圖1-4分別是4種不同模板信號：方波信號、正弦信號、指數(shù)衰減信號、指數(shù)衰減正弦信號的的數(shù)據(jù)圖、信號圖和線性相關圖； 比較四個不

8、同模板信號圖可以看出：四種信號的噪聲信號圖在相同區(qū)間0 1000 -5 5變化相似；方波和指數(shù)衰減的有用信號圖在相同區(qū)間0 1000 0 4變化相似，峰值的范圍和取值相似，正弦和指數(shù)衰減正弦信號的有用信號圖在相同區(qū)間0 1000 -5 5變化相似，方波信號、正弦信號、指數(shù)衰減正弦信號的觀測信號圖在相同區(qū)間0 1000 -10 10變化相似，指數(shù)衰減信號在區(qū)間0 1000 -5 5變化；方波、正弦信號的線性相關圖在相同區(qū)間變化；指數(shù)衰減信號、指數(shù)衰減正弦信號的變化趨勢相似，但是峰值等不相同，四種信號的橫坐標取值都是-1000 1000，縱坐標各不相同，總體的圖形變化趨勢是類似的。思考題：嘗試修改

9、程序，包括改變仿真信號中模版的形狀，噪聲的強弱，噪聲的類型（對白噪聲濾波可以獲得各種有色噪聲），哪些因素會影響相關函數(shù)的結果？clear;clc; %clear清空內存變量,clc只清空工作區(qū)內的內容np = 0:99; %從0到99取值%load ecgdata %導入心電信號%load eegdata %導入腦電信號%load icpdata %導入顱內壓信號%load respdata %導入呼吸信號%p=(ecgdata(1:100); %產(chǎn)生100個心電信號的數(shù)據(jù)% p=(eegdata(1:100); %產(chǎn)生100個腦電信號的數(shù)據(jù)% p=(icpdata(1:100); %產(chǎn)生10

10、0個顱內壓信號的數(shù)據(jù)% p=(respdata(1:100); %產(chǎn)生100個呼吸信號的數(shù)據(jù)%p = sin(pi/5*np); % 正弦% p = exp(-0.06*np); % 指數(shù)衰減% p = sin(pi/5*np).*exp(-0.06*np); % 指數(shù)衰減正弦 p = ones(size(np); % 方波figure; %畫圖subplot(1,1,1); plot(np,p); %數(shù)據(jù)圖n1=1:100load ecgdata; %鏈接心電信號 x1 = ecgdata(n1); %產(chǎn)生100個心電信號的數(shù)據(jù)load eegdata; %鏈接腦電信號 x2 = eegda

11、ta (n1); %產(chǎn)生100個腦電信號的數(shù)據(jù)load icpdata; %鏈接顱內壓信號 x3 = icpdata (n1); %產(chǎn)生100個顱內壓信號的數(shù)據(jù)load respdata; %鏈接呼吸信號 x4 = respdata (n1); %產(chǎn)生100個呼吸信號的數(shù)據(jù)figure;subplot(2,2,1); plot(n1,x1);title(心電信號波形); %畫出100點長的心電信號波形圖subplot(2,2,2); plot(n1,x2);title(腦電信號波形); %畫出100點長的腦電信號波形圖subplot(2,2,3); plot(n1,x3);title(顱內壓信

12、號波形);%畫出100點長的顱內壓信號波形圖subplot(2,2,4); plot(n1,x4);title(呼吸信號波形); %畫出100點長的呼吸信號波形圖n = 1:1000; %定義1000點長的隨機信號w = randn(size(n); %產(chǎn)生1000點長的隨機信號噪聲%w = 0.5*randn(size(n); %產(chǎn)生1000點長的均值為0，方差為0.5隨機信號噪聲%w = randn(size(n); %產(chǎn)生1000點長的均值為0，方差為1隨機信號噪聲%w = 2*randn(size(n); %產(chǎn)生1000點長的均值為0，方差為2隨機信號噪聲s = zeros(size(

13、n); %產(chǎn)生1000點長的有用信號A = 3; % 衰減系數(shù)%A = 1; % 衰減系數(shù)為1%A = 6; % 衰減系數(shù)為6%A = 9; % 衰減系數(shù)為9s(100:199) = s(100:199)+A*p;s(500:599) = s(500:599)+A/3*p;s(800:899) = s(800:899)+A/3/3*p;x = s+w; %仿真回波與白噪聲疊加構造仿真信號figure; %畫圖subplot(3,1,1); plot(n,w); title(Noise); %作出噪聲信號的圖，并添加標題subplot(3,1,2); plot(n,s); title(Signa

14、l); %作出有用信號的圖，并添加標題subplot(3,1,3); plot(n,x); title(Signal with Noise); %作出觀測信號圖，并添加標題y=medfilt1(w,100); %用中值濾波函數(shù)medfilt1對噪聲的100個點取中值figure; %畫圖subplot(3,1,1);plot(n,w);legend(orignal signal); %建立原始信號的子圖 subplot(3,1,2);plot(n,y);title(N=60);legend(中值濾波);%建立子圖，添加標題y1=imfilter(w,100);subplot(3,1,3);pl

15、ot(n,y1);title(N=60);legend(均值濾波);%建立子圖，添加標題p = p,zeros(1,length(x)-length(p); % 如果要求歸一化相關系數(shù)（相干系數(shù)），兩個序列要同樣長Rps = xcorr(s,p,coeff);%對信號做自相關函數(shù)Rpw = xcorr(w,p,coeff);%對噪聲信號做自相關函數(shù)Rpx = xcorr(x,p,coeff);%對觀測信號做自相關函數(shù)%線性相關n2 = (n(1)-n(end):(n(end)-n(1); %線性相關的范圍figure; %畫圖subplot(3,1,1); plot(n2,Rpw); titl

16、e(Rpw of p(n) and w(n);subplot(3,1,2); plot(n2,Rps); title(Rps of p(n) and s(n); subplot(3,1,3); plot(n2,Rpx); title(Rpx of p(n) and x(n);%線性卷積Cpw=conv(w,p);Cps=conv(s,p);Cpx=conv(x,p);figure;subplot(3,1,1); plot(Cpw);xlim(0,2001); title(Cpw of p(n) and w(n);subplot(3,1,2); plot(Cps);xlim(0,2001); t

17、itle(Cps of p(n) and s(n); subplot(3,1,3); plot(Cpx);xlim(0,2001); title(Cpx of p(n) and x(n);%線性相干Rps2 = xcorr(s,p,coeff);Rpw2 = xcorr(w,p,coeff);Rpx2 = xcorr(x,p,coeff);figure;subplot(3,1,1); plot(Rpw2); axis(0,2001,-1,1); title(Rpw2 of p(n) and w(n);subplot(3,1,2); plot(Rps2); axis(0,2001,-1,1);

18、title(Rps2 of p(n) and s(n); subplot(3,1,3); plot(Rpx2); axis(0,2001,-1,1); title(Rpx2 of p(n) and x(n);%循環(huán)相關Vpw=circlel(p);Rpw1=w*Vpw;Vps=circlel(p);Rps1=s*Vps;Vpx=circlel(p);Rpx1=x*Vpx;figure;subplot(3,1,1); plot(Rpw1);xlim(0,1000); title(Rpw1 of p(n) and w(n);subplot(3,1,2); plot(Rps1);xlim(0,100

19、0); title(Rps1 of p(n) and s(n); subplot(3,1,3); plot(Rpx1);xlim(0,1000); title(Rpx1 of p(n) and x(n);%循環(huán)卷積Vpw=circler(p);Cpw1=w*Vpw;Vps=circler(p);Cps1=s*Vps;Vpx=circler(p);Cpx1=x*Vpx;figure;subplot(3,1,1); plot(Cpw1);xlim(0,1000); title(Cpw1 of p(n) and w(n);subplot(3,1,2); plot(Cps1);xlim(0,1000)

20、; title(Cps1 of p(n) and s(n); subplot(3,1,3); plot(Cpx1);xlim(0,1000); title(Cpx1 of p(n) and x(n); 圖5 方波 圖6 四種信號波形圖 圖7 方波的噪聲、有用信號和觀測信號圖 圖8 對噪聲的中值濾波和均值濾波 圖9 方波的線性相關圖 圖10 方波的線性卷積圖 圖11 方波的線性相干圖 圖12 方波的循壞相關圖圖13 方波的循壞卷積圖答：只改變仿真信號中模版的形狀：通過對p的注釋與選擇，分別選擇8種不同模板信號：方波信號、正弦信號、指數(shù)衰減信號、指數(shù)衰減正弦信號、心電信號、腦電信號、顱內壓信號、呼

21、吸信號；A=3，噪聲均值為0，方差為1.np = 0:99; %從0到99取值%load ecgdata %導入心電信號%load eegdata %導入腦電信號%load icpdata %導入顱內壓信號%load respdata %導入呼吸信號%p=(ecgdata(1:100); %產(chǎn)生100個心電信號的數(shù)據(jù)% p=(eegdata(1:100); %產(chǎn)生100個腦電信號的數(shù)據(jù)% p=(icpdata(1:100); %產(chǎn)生100個顱內壓信號的數(shù)據(jù)% p=(respdata(1:100); %產(chǎn)生100個呼吸信號的數(shù)據(jù)%p = sin(pi/5*np); % 正弦% p = exp(-

22、0.06*np); % 指數(shù)衰減% p = sin(pi/5*np).*exp(-0.06*np); % 指數(shù)衰減正弦p = ones(size(np); % 方波方波的線性相關和卷積相關圖14方波的線性相關和卷積相關圖正弦波的線性相關和卷積相關 圖15正弦波的線性相關和卷積相關圖 指數(shù)衰減波的線性相關和卷積相關 圖16指數(shù)衰減波的線性相關和卷積相關圖 指數(shù)正弦衰減波的線性相關和卷積相關 圖17 指數(shù)正弦衰減波的線性相關和卷積相關圖心電波的線性相關和卷積相關 圖18心電波的線性相關和卷積相關圖 腦電波的線性相關和卷積相關 圖19腦電波的線性相關和卷積相關圖 顱內壓波的線性相關和卷積相關 圖20

23、顱內壓波的線性相關和卷積相關圖 呼吸波的線性相關和卷積相關 圖21呼吸波的線性相關和卷積相關圖答：通過圖像觀察，只改變仿真信號中模版的形狀，產(chǎn)生的相關函數(shù)圖的幅值不同，所以它會影響線性相關函數(shù)的結果，其中線性相干是線性相關的歸一化，xym=rxymrxx0ryy01/2形狀相同，故略去，它們的結論也相同；相關的取值范圍是-1000 1000,卷積的范圍是0 2000.只改變噪聲的強弱：A=3，分別選擇均值都為0，方差為0.5,1,2；方波的不同噪聲強度的線性相關圖和卷積相關圖%w = 0.5*randn(size(n); %產(chǎn)生1000點長的均值為0，方差為0.5隨機信號噪聲 圖22 均值都為

24、0，方差為0.5的方波線性相關圖和卷積相關圖%w = randn(size(n); %產(chǎn)生1000點長的均值為0，方差為1隨機信號噪聲 圖23 均值都為0，方差為1的方波線性相關圖和卷積相關圖%w = 2*randn(size(n); %產(chǎn)生1000點長的均值為0，方差為2隨機信號噪聲 圖23 均值都為0，方差為2的方波線性相關圖和卷積相關圖答：從圖看出，方差為0.5,1,2，噪聲強度越大，方波線性相關的峰值越小，卷積相關的峰值變化不大，會影響結果；正弦波不同噪聲強度的線性相關圖和卷積相關圖A=3，分別選擇均值都為0，方差為0.5,1,2；正弦波的均值為0，方差為0.5隨機信號噪聲 圖24均值

25、都為0，方差為0.5的正弦波線性相關圖和卷積相關圖正弦波的均值為0，方差為1隨機信號噪聲 圖25均值都為0，方差為1的正弦波線性相關圖和卷積相關圖正弦波的均值為0，方差為2隨機信號噪聲 圖26均值都為0，方差為2的正弦波線性相關圖和卷積相關圖答：從圖看出，方差為0.5,1,2，噪聲強度越大，正弦波線性相關的峰值越小，卷積相關峰值變化不大，會影響結果；心電波的不同噪聲強度的線性相關圖和卷積相關圖A=3，分別選擇均值都為0，方差為0.5,1,2；心電波的均值為0，方差為0.5隨機信號噪聲 圖27 均值都為0，方差為0.5的心電波線性相關圖和卷積相關圖心電波的均值為0，方差為1隨機信號噪聲 圖28

26、均值都為0，方差為1的心電波線性相關圖和卷積相關圖心電波的均值為0，方差為2隨機信號噪聲 圖29 均值都為0，方差為2的心電波線性相關圖和卷積相關圖答：從圖看出，方差為0.5,1,2，噪聲強度越大，心電線性相關的峰值越小，卷積相關的峰值變化不大，會影響結果；只改變噪聲的類型（對白噪聲濾波可以獲得各種有色噪聲）：對噪聲分別進行中值濾波和均值濾波，A=3，均值為0，方差為1未濾波的方波的線性相關和卷積相關 圖30 未濾波的方波的線性相關和卷積相關圖y=medfilt1(w,100); %中值濾波 下圖為中值濾波的方波的線性相關和卷積相關 圖31中值濾波的方波的線性相關和卷積相關圖y=imfilte

27、r(w,100); %均值濾波 下圖為均值濾波的方波的線性相關和卷積相關 圖32 均值濾波的方波的線性相關和卷積相關圖答：從圖看出，改變噪聲的類型，相關函數(shù)的變化趨勢相似，所以不影響相關函數(shù)的結果；只改變衰減系數(shù)A為：1,3,9，模板為方波，均值為0，方差為1的2種線性函數(shù)圖A = 1; % 衰減系數(shù)為1 圖33 衰減系數(shù)為1的方波線性相關和卷積相關圖A = 3; % 衰減系數(shù)為3 圖34 衰減系數(shù)為3的方波線性相關和卷積相關圖A = 9; % 衰減系數(shù)為9 圖35 衰減系數(shù)為9的方波線性相關和卷積相關圖答：從圖中可以看出，隨著衰減系數(shù)1,3,9的增大，線性相關和卷積相關函數(shù)圖的峰值越大；卷積

28、之前先對信號進行function調用；%circlel循壞相關函數(shù) %circler循壞卷積函數(shù)function v=circlel(y) function v=circler(y)N=length(y); N=length(y);v=zeros(N,N); v=zeros(N,N);for i=1:N for i=1:N for j=1:N for j=1:N v(i,j)=y(j); v(i,j)=y(j); end endL=y(1); L=y(N);for k=1:N-1 for k=N:-1:2 y(k)=y(k+1); y(k)=y(k-1);end endy(N)=L; y(1)=L;end end v=v;模板信號為方波，A=3，噪聲均值為0，方差為1的2種循環(huán)