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

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)=L;end end v=v;模板信号为方波，A=3，噪声均值为0，方差为1的2种循环