数字相关和数字卷积 实验报告.doc

《数字相关和数字卷积 实验报告.doc》由会员分享，可在线阅读，更多相关《数字相关和数字卷积 实验报告.doc（25页珍藏版）》请在一课资料网上搜索。

1、生物医学信号处理评分大理大学实验报告 课程名称 生物医学信号处理 实验名称 数字相关和数字卷积 专业班级 姓 名 羽卒兰cl 学 号 实验日期 实验地点 20152016学年度第 3 学期一、 实验目的熟悉数字相关的运算，初步在信号处理中应用相关技术。二、实验环境 1、硬件配置：Intel（R） Core(TM) i5-4210U CPU 1.7GHz 1.7GHz 安装内存（RAM）：4.00GB 系统类型：64位操作系统 2、软件环境：MATLAB R2013b软件三、实验内容(包括本实验要完成的实验问题及需要的相关知识简单概述)已知发射波形，利用相关技术，在有强背景噪声的情况下检测回波的

2、延时和强度。首先使用已知信号模版及其若干次衰减延迟生成仿真回波波形，然后与白噪声背景叠加，构造仿真信号。然后计算模版与仿真信号的相关函数，判断回波位置及相对强度。四、实验结果与分析 （包括实验原理、数据的准备、运行过程分析、源程序（代码）、图形图象界面等）实验原理相关可以从时域角度表现信号间的相似（关联）程度，是统计信号处理最基本的手段之一。设有离散信号x(n)和y(n)，线性相关函数定义为 实际采集的信号总是有限长度，用有限的样本估计相关（自相关）函数 求和项总数不是N而是N-|m|，因为当n=N-|m|-1时，n+|m|=N-1。此时xn+m已经到了数据边沿。这种估计是渐进无偏估计和一致估

3、计。计算中，只要将其中一个序列反转，就可以用计算线性卷积的程序计算线性相关 因此可以用FFT来加速相关运算，即对序列补零后，用循环相关计算线形相关，然后用循环卷积的快速算法计算循环相关，得到最终结果。 源代码：% 估计两个相似信号间的时间延迟np = 0:99; %从0到99取值% p = sin(pi/5*np); % 正弦% p = exp(-0.06*np); % 指数衰减% p = sin(pi/5*np).*exp(-0.06*np); % 指数衰减正弦p = ones(size(np); % 方波figure;subplot(2,2,1); plot(np,p);%数据图n = 1

4、1000; %定义1000点长的随机信号w = randn(size(n); %产生1000点长的随机信号噪声s = zeros(size(n); %产生1000点长的有用信号%通过改变衰减系数对已知信号模版进行若干次衰减，生成仿真回波A = 3; % 衰减系数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); %如果要求归一化相关系数（相干系数），两个序列要同样长Rps = xcorr(s,p,coeff); %对有用信号做互相关函数Rpw = xcorr(w,p,coeff); %对噪声信号做互相关函数Rpx = xcorr(x,p,coeff); %对观

6、测信号做互相关函数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 方波信号的数据图a、信号图b和线性相关图c2. 模板为正弦信号，A=3，噪声均值为0，方差为1a） b）c） 图2 正弦信号的数据图a、信号图b和线性相关图c3. 模板为指数衰减信号，A=3，噪声均值为0，方差为1a） b） c） 图3 指数衰减信号的数据图a、信号图b和线性相关图c4. 模板为指数衰减正弦信号，A=3，噪声均值为0，方差为1a） b） c） 图4 指数衰减正弦信号的数据图a、信号图b和线性相关图c分析： 图1-4分别是4种不同模板信号：方波信号、正弦信号、指数衰减信号、指数衰减正弦信号的的数据图、信号图和线性相关图； 比较四个不

8、同模板信号图可以看出：四种信号的噪声信号图在相同区间0 1000 -5 5变化相似；方波和指数衰减的有用信号图在相同区间0 1000 0 4变化相似，峰值的范围和取值相似，正弦和指数衰减正弦信号的有用信号图在相同区间0 1000 -5 5变化相似，方波信号、正弦信号、指数衰减正弦信号的观测信号图在相同区间0 1000 -10 10变化相似，指数衰减信号在区间0 1000 -5 5变化；方波、正弦信号的线性相关图在相同区间变化；指数衰减信号、指数衰减正弦信号的变化趋势相似，但是峰值等不相同，四种信号的横坐标取值都是-1000 1000，纵坐标各不相同，总体的图形变化趋势是类似的。思考题：尝试修改

9、程序，包括改变仿真信号中模版的形状，噪声的强弱，噪声的类型（对白噪声滤波可以获得各种有色噪声），哪些因素会影响相关函数的结果？clear;clc; %clear清空内存变量,clc只清空工作区内的内容np = 0:99; %从0到99取值%load ecgdata %导入心电信号%load eegdata %导入脑电信号%load icpdata %导入颅内压信号%load respdata %导入呼吸信号%p=(ecgdata(1:100); %产生100个心电信号的数据% p=(eegdata(1:100); %产生100个脑电信号的数据% p=(icpdata(1:100); %产生10

10、0个颅内压信号的数据% p=(respdata(1:100); %产生100个呼吸信号的数据%p = sin(pi/5*np); % 正弦% p = exp(-0.06*np); % 指数衰减% p = sin(pi/5*np).*exp(-0.06*np); % 指数衰减正弦 p = ones(size(np); % 方波figure; %画图subplot(1,1,1); plot(np,p); %数据图n1=1:100load ecgdata; %链接心电信号 x1 = ecgdata(n1); %产生100个心电信号的数据load eegdata; %链接脑电信号 x2 = eegda

11、ta (n1); %产生100个脑电信号的数据load icpdata; %链接颅内压信号 x3 = icpdata (n1); %产生100个颅内压信号的数据load respdata; %链接呼吸信号 x4 = respdata (n1); %产生100个呼吸信号的数据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); %产生1000点长的随机信号噪声%w = 0.5*randn(size(n); %产生1000点长的均值为0，方差为0.5随机信号噪声%w = randn(size(n); %产生1000点长的均值为0，方差为1随机信号噪声%w = 2*randn(size(n); %产生1000点长的均值为0，方差为2随机信号噪声s = zeros(size(

13、n); %产生1000点长的有用信号A = 3; % 衰减系数%A = 1; % 衰减系数为1%A = 6; % 衰减系数为6%A = 9; % 衰减系数为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); %用中值滤波函数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); % 如果要求归一化相关系数（相干系数），两个序列要同样长Rps = xcorr(s,p,coeff);%对信号做自相关函数Rpw = xcorr(w,p,coeff);%对噪声信号做自相关函数Rpx = xcorr(x,p,coeff);%对观测信号做自相关函数%线性相关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);%循环相关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);%循环卷积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种不同模板信号：方波信号、正弦信号、指数衰减信号、指数衰减正弦信号、心电信号、脑电信号、颅内压信号、呼

21、吸信号；A=3，噪声均值为0，方差为1.np = 0:99; %从0到99取值%load ecgdata %导入心电信号%load eegdata %导入脑电信号%load icpdata %导入颅内压信号%load respdata %导入呼吸信号%p=(ecgdata(1:100); %产生100个心电信号的数据% p=(eegdata(1:100); %产生100个脑电信号的数据% p=(icpdata(1:100); %产生100个颅内压信号的数据% p=(respdata(1:100); %产生100个呼吸信号的数据%p = sin(pi/5*np); % 正弦% p = exp(-

22、0.06*np); % 指数衰减% p = sin(pi/5*np).*exp(-0.06*np); % 指数衰减正弦p = ones(size(np); % 方波方波的线性相关和卷积相关图14方波的线性相关和卷积相关图正弦波的线性相关和卷积相关 图15正弦波的线性相关和卷积相关图 指数衰减波的线性相关和卷积相关 图16指数衰减波的线性相关和卷积相关图 指数正弦衰减波的线性相关和卷积相关 图17 指数正弦衰减波的线性相关和卷积相关图心电波的线性相关和卷积相关 图18心电波的线性相关和卷积相关图 脑电波的线性相关和卷积相关 图19脑电波的线性相关和卷积相关图 颅内压波的线性相关和卷积相关 图20

23、颅内压波的线性相关和卷积相关图 呼吸波的线性相关和卷积相关 图21呼吸波的线性相关和卷积相关图答：通过图像观察，只改变仿真信号中模版的形状，产生的相关函数图的幅值不同，所以它会影响线性相关函数的结果，其中线性相干是线性相关的归一化，xym=rxymrxx0ryy01/2形状相同，故略去，它们的结论也相同；相关的取值范围是-1000 1000,卷积的范围是0 2000.只改变噪声的强弱：A=3，分别选择均值都为0，方差为0.5,1,2；方波的不同噪声强度的线性相关图和卷积相关图%w = 0.5*randn(size(n); %产生1000点长的均值为0，方差为0.5随机信号噪声 图22 均值都为

24、0，方差为0.5的方波线性相关图和卷积相关图%w = randn(size(n); %产生1000点长的均值为0，方差为1随机信号噪声 图23 均值都为0，方差为1的方波线性相关图和卷积相关图%w = 2*randn(size(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 均值滤波的方波的线性相关和卷积相关图答：从图看出，改变噪声的类型，相关函数的变化趋势相似，所以不影响相关函数的结果；只改变衰减系数A为：1,3,9，模板为方波，均值为0，方差为1的2种线性函数图A = 1; % 衰减系数为1 图33 衰减系数为1的方波线性相关和卷积相关图A = 3; % 衰减系数为3 图34 衰减系数为3的方波线性相关和卷积相关图A = 9; % 衰减系数为9 图35 衰减系数为9的方波线性相关和卷积相关图答：从图中可以看出，随着衰减系数1,3,9的增大，线性相关和卷积相关函数图的峰值越大；卷积

28、之前先对信号进行function调用；%circlel循坏相关函数 %circler循坏卷积函数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)

29、L;end end v=v;模板信号为方波，A=3，噪声均值为0，方差为1的2种循环函数： 图36 1001点循坏相关函数 图37 1001点循坏卷积函数Vpw=circlel(p);Rpw=w*Vpw;Vps=circlel(p);Rps=s*Vps;Vpx=circlel(p);Rpx=x*Vpx;Vpw=circler(p);Cpw=w*Vpw;Vps=circler(p);Cps=s*Vps;Vpx=circler(p);Cpx=x*Vpx;模板信号为正弦波，A=3，噪声均值为0，方差为1的2种循环函数： 图38 1001点循坏相关函数 图39 1001点循坏卷积函数Vpw=circ

30、ler(p);Cpw=w*Vpw;Vps=circler(p);Cps=s*Vps;Vpx=circler(p);Cpx=x*Vpx;Vpw=circlel(p);Rpw=w*Vpw;Vps=circlel(p);Rps=s*Vps;Vpx=circlel(p);Rpx=x*Vpx;答：循坏相关就是周期序列向左移动一个位置，循坏卷积就是周期序列向右移动一个位置。改变模板信号，噪声的循环函数变化相似，有用信号和观测信号的循环函数变化改变；只改变噪声的强弱：A=3，模板为方波，分别选择均值都为0，方差为0.5,1,2；方波的均值为0，方差为0.5随机信号噪声 图40均值为0，方差为0.5的方波的循

31、环相关和循环卷积图方波的均值为0，方差为1随机信号噪声 图41均值为0，方差为1的方波的循环相关和循环卷积图方波的均值为0，方差为2随机信号噪声 图42均值为0，方差为2的方波的循环相关和循环卷积图答：从图中可以看出，噪声的方差分别为0.5,1,2，噪声强度越大，循环函数的峰值越大。五、实验小结： （包括主要实验问题的最终结果描述、详细的收获体会，待解决的问题等） 卷积之前先对信号进行function调用；调用的function函数，导入的各种信号等要选择在同一文件目录下进行运行，不然会出错；进行卷积等运算的序列的长度要一致，不然会报错；S和W的维数要相同，都是1*1001，不然无法进行运算；要确定各种函数图的取值范围，在同一范围内进行比较；判断各因素的影响的时候要控制单一变量，只改变一种因素，其它因素不变才能进行判断。 手写签名： 第25页 共25页