医学信号处理实验报告-两路信号关系衡量

电子科技大学生命科学与技术学院标准实验报告
（实验）课程名称生物医学信号处理
2018-2019-第2学期
电子科技大学教务处制表
学生姓名学号
指导教师实验时间2019.4.16
一、实验室名称：品学楼B302
二、实验名称：两路信号间的关系衡量
三、实验学时：2
四、实验原理：
1.信噪比(signal-to-noise ratio):是描述信号中有效成分与噪声成分的比例关系参数，单位为dB。

假设不含噪声的信号为，外加噪声以后的信号为，则信号的信噪比定义为
其中代表信号的方差
在给定信噪比的情形下，要求解系数，则其计算公式为
2.皮尔逊相关系数
在统计学中，皮尔逊相关系数(Pearson correlation coefficient)，通常用R或表示，是用来度量两个变量X和Y之间的相互关系(线性相关)的，取值范围在 [-1,+1] 之间。

它在学术研究中被广泛应用来度量两个变量线性相关性的强弱。

在作为衡量线性回归效果时，常使用
对于随机变量X和Y，皮尔森相关系数的求解公式为：
其中代表X与Y的协方差，和分别代表X 和Y的方差。

当相关性为 1 时，X与Y的关系可以表示为，其中a>0；当相关性为 -1 时，X与Y的关系可以表示为，其中a<0。

如果X与Y相互独立，那么相关性为0。

两个变量为正相关，则皮尔逊相关系数在0与1之间，两个变量为负相关，则皮尔逊相关系数在-1与0之间。

3.自相关检测含噪信号周期:
若为周期性的随机信号，为随机噪声信号，为实际接收到的信号，则。

因此
因此接收到的信号自相关可以分解为四个与发送信号和噪声有关的自相关函数，其中，的自相关函数只在零点处有最大值，其余点可大致认为0。

而由于信号和噪声的独立性，
而噪声信号的期望值为0，所以在零点之外，=0，因而，由于是周期信号
所以有
因此自相关函数的周期与原始信号一致。

4.延时相关：
5.频域相干
设有两个信号，它们的频域相干函数(幅值平方相干函数)定义如下
式中，表示两个信号之间的互功率谱，即这两个信号的循环互相关函数的离散傅里叶变换，
为各自的功率谱，这里的代表频率， 的取值范围为0~1之间。

当它等于1时，说明两个信号是完全相干的，即一个信号可以完全由另一个信号决定；当它等于0时，说明这两个信号完全不相干，即这两个信号完全独立；当它在0~1之间时，说明这两个信号存在部分相干性，即非线性关系或有外界的干扰存在。

6.线性相关 设有离散信号，则其线性相关函数为
[]00(0)(1)(1)0(0)(1)(2)0(0)(1)(1)(1)(2)00
(1)(0)(2)(1)00y y y N y y y x x x N y y y N y y N y N -⎡⎤⎢⎥⎢⎥⎢⎥=-⎢⎥-⎢⎥⎢⎥--⎣⎦
7.Mscohere 函数
1.使用mscohere 函数的时候为了使图像更平滑，可以使用更短的窗（其实也代表着分段周期图法分段长度更短），主瓣更宽窗函数。

2.Noverlap的作用
因为加窗丢失了两端的信息，所以分段时相邻段之间要有重复的部分，一般为1/2或1/3的窗长。

效果如下图：
五、实验目的：
理解相关系数，序列相关，频率相干的方法，应用：
1.了解噪声对信号检测的影响。

2.了解利用延时相关估计距离的方法。

3.利用频率相干分析波形频率成分。

4.了解线性相关的计算。

六、实验内容：
上机题1：成对相关
1、在波形产生函数中选sinc波形发生器（s）和randn白噪声发生器（n），叠加成一个信号（x=s+weight*n），请改变白噪声信号的加权值，使得信噪比约为10,5,0，-5，-10 dB，求得5个weight 值。

画出s和x信号（一张图），观察。

2、用s信号与五个x信号，求他们的成对相关，corr或用corrcoef 函数，画出相关系数与信噪比的关系图，文字描述图说明了什么。

3、用s信号与x=weight*n-s做成对相关，画出相关系数与信噪比的关系图，观察s与x信号，文字描述图说明了什么。

4、结合上述结果，请问相关系数为正和负分别代表了什么含义？相关系数值大小代表了什么含义？
上机题2：相关——信号的检测
1、任意产生一个白噪信号(w) 和一个双频率的信号(s) (for example, s = 2sin(w1n)+4sin(w2n))
2、信号s的周期是多少？仿真 s, w, and x, (x=s+a*w，a=1) , 计算三个信号的线性相关，并画出时间信号和线性相关结果。

3、你能否从线性相关函数中检测出有用信号的周期？报导周期值。

4、改变a的大小，报导x与s的相关系数的变化情况，也即改变信噪比（信号能量除以噪声能量），相似程度变化情况图。

上机题3：延时相关——估计距离的一种方法
两个拾音器接收到的声波信号，分别来自一个独立声源的两条路径：直达声和反射声。

声传播的速度为，衰减系数参数为，系统的两个输出分别描述为：
X(t) = p(t)+ A1*p(t -)
y(t) = A*p(t - )+ A2*p(t - )
预测：R xy()在、出现极大值如下图：
上机题4：频域相干——脑电信号的频域相干
1、学习使用mscohere函数，估计两个信号频域的关系。

2、任意产生两个信号，画出两个信号分别的功率谱，及他们的频域相干图，观察相干大小和两个信号功率谱之间的关系，请描述。

3、用睁眼和闭眼脑电信号（数据文件eegclose.mat 和eegopen.mat ，Fs=250 Hz，幅度单位：微伏），用自己信号（学号
*班号）与其他任意一路信号做频域相干，观察睁眼和闭眼的频域相干的差异（4个频段，列表方式），并描述。

4、计算自己信号的闭眼和睁眼信号的频域相干（4个频段，列表方式），描述结果。

5、总结睁眼闭眼脑电信号的频域关系？
上机题5：编制函数实现两个随机序列的线性相关与MATLAB 自带函数xcorr 结果进行比较。

10()()()
N xy n r m x n y n m -==+∑
七、实验器材（设备、元器件）： matlab2014b
八、实验步骤：构思，写代码，写报告，修改代码，完成报告。

九、实验数据及结果分析：
（一）程序：
t=linspace(-5,5,1000);
s = sinc(t);
n = randn(1,length(t));
snr = linspace(10,-10,5);
snr2 = linspace(50,-50,101);
snr_len = length(snr);
snr_len2 = length(snr2);
snrt = zeros(1,snr_len);
weight = zeros(1,snr_len);
x = zeros(1000,snr_len);
y = zeros(1000,snr_len2);
corr_sx = zeros(snr_len2,1);
for i = 1:length(snr)
weight(i) = sqrt(var(s)/(var(n)*10^(snr(i)/10))); x(:,i) = s + weight(i)*n ;
figure(1)
subplot(5,1,i)
plot(x(:,i)); hold on;plot(s);hold on;
title(['信噪比为',num2str(snr(i))])
snrt(i) = 10*log10(var(s)/var(n)/(weight(i)^2)); corr_sx(i) = corr(s',x(:,i));
hold on
end
weight
for i = 1:length(snr2)
weight(i) = sqrt(var(s)/(var(n)*10^(snr2(i)/10))); y(:,i) = weight(i)*n-s ;
snrt(i) = 10*log10(var(s)/var(n)/(weight(i)^2)); corr_sx(i) = corr(s',y(:,i));
hold on
end
figure(2)
plot(snr2,corr_sx)
xlabel('信噪比/dB')
ylabel('相关系数')
结果：
Weight：0.0906 0.1610 0.2864 0.5093 0.9056（对应10 ，5 ，0 ，-5 ，-10）
图1 不同信噪比的波形
图2 相关系数随信噪比变化图1
图3 相关系数随信噪比变化图2 (二)
程序：
clear all
clc
t=linspace(-5,5,1000);
s = 2*sin(2*pi*t)+4*sin(pi*t);
T = 2*pi/pi*100;
n = randn(1,length(t));
a = linspace(0,30,31);
a_len = length(a);
a_t = zeros(1,a_len);
weight = zeros(1,a_len);
x = zeros(1000,a_len);
corr_sx = zeros(a_len,1);
xcorr_x = zeros(2*1000-1,a_len);
for i = 1:length(a)
x(:,i) = s + a(i)*n ;
corr_sx(i) = corr(s',x(:,i));
if a(i)==1
figure(1)
plot(x(:,i),'y');
hold on;
plot(s,'r');
hold on;
figure(2)
xcorr_x(:,i) = xcorr(x(:,i),s,'biased')'; [pks,locs]=findpeaks(xcorr_x(:,i));
lo_cs = find(pks>0.1);
locs = locs(lo_cs);
T_c = mean(diff(locs));
plot( xcorr_x(:,i));
end
end
figure(3)
plot(0:30,corr_sx,'o');
disp('实际周期为：');disp(T);
disp('相关测出周期为：');disp(T_c);
结果：
实际周期为：200 相关测出周期为：199.8750
图4 双频信号与含有噪声的双频信号
图5 含有/不含有噪声的信号线性相关结果（有偏）
图6 相关系数与系数a关系图
（三）
程序：
clear all;
N=1000; %长度
Fs=50; %采样频率
n=0:N-1; t=n/Fs; %时间序列
A=0.4;A1=0.5;A2=0.6; %衰减系数
c0=340; %c0
d1=620
d2=500 %请自己给两个距离的参数，不要太小即可t1=d1/c0;
t2=(d1+2*d2)/c0;
tc=2*(d1+d2)/c0;
Lag=500; %最大延迟样点数
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));
figure(2)
plot((1:N)/Fs,yt);
[Rpp,lags] = xcorr(pt,'biased'); %p(t)自相关
[Rxx,lagx] = xcorr(xt,'biased'); %x(t)自相关
[Ryy,lagt] = xcorr(yt,'biased'); %y(t)自相关
[Rxy,lagy] = xcorr(xt,yt,'biased'); %x(t)与y(t)互相关
rts=lags/Fs;rtx=lagx/Fs;rtt=lagt/Fs;rty=lagy/Fs;
figure(3)
[pks,locs]=findpeaks(Rxy);
lo_cs = find(pks>1E-3);
locs = (locs(lo_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)自相关')
xlabel('时间/s')
subplot(2,2,2)
plot(rtx,Rxx)
title('x(t)自相关')
xlabel('时间/s')
subplot(2,2,3)
plot(rtt,Ryy)
title('%y(t)自相关')
xlabel('时间/s')
subplot(2,2,4)
plot(rty,Rxy)
title('%x(t)与y(t)互相关')
xlabel('时间/s')
d1_c = (locs(3) * c0 - locs(2) * c0)/2
d2_c = ((locs(4) * c0 - d1_c)/2 + (-locs(1) * c0 - d1_c)/2)/2结果：
图7 三个信号自相关及x，y互相关图
d1d2d1（测）d2（测）
620500612499.8
650420646418.2
800700792.2702.1
420650411.4651.1
表1 d1，d2测量情况表
（四）
程序：
clear all ;
close=load('eegclose');
open=load('eegopen');
close1 = close.eegclose(:,8*2);
close2 = 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),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,'XTickmode','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',locate_mid)
结果：
图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 同路信号相干结果表
（五）
程序：
clear all;
clc
a = [1 2 5 6 7];
b = [2 5 6];
s1 = Icorr(a,b)
s2 = xcorr(a,b)
function [s] = Icorr(s1,s2)
% 计算两个函数的线性相关
s1_s = size(s1);
s2_s = size(s2);
if s1_s(1) ~= 1
s1 = s1';
end
if s2_s(1) ~= 1
s2 = s2';
end
s1_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,zeros(1,dlen)];
end
y = zeros(s_len,s_len*2-1);
for i = 1:s_len
y(i,s_len+i-1:-1:i) = s1;
end
s = s2*y;
s = s(1:s_len*2-1-dlen);
结果：
Icorr： 14 47 82 65 42 17 6
Xcorr：0.0000 0.0000 6.0000 17.0000 42.0000 65.0000 82.0000 47.0000 14.0000
十、总结及心得体会：
（一）
1. weight值随着信噪比减小而升高。

2.构造信号（x=s+weight*n），信噪比越高相关系数越大。

在20dB 时已经接近最大为1，在-20dB时接近最小值0.04。

构造信号（x=weight*n-s）时，在-20dB时已经接近最小为-1，在20dB时接近最大值-0.04。

这说明，如果要从含噪声的信号中检测到信号，一般手段的极限在[-20,20]dB之间。

3.相关系数的符号代表了进行相关计算两信号序列的变化趋势关系。

如为正，说明一个信号序列中增大时，另一个相应的增大，反之减小，称为具有正相关性。

符号为负则相反，称为具有负相关性。

（二）
1.将信号与含噪声的信号进行线性相关可以有效地识别出周期。

设想该性质是由于加入的是randn函数生成的白噪声信号。

改为加入rand函数生成的噪声，自相关的方法仍能准确体现出信号的周期性。

变化的只有相关系数的变化趋势变得趋于线性，不过这一点在意料之中。

2.如果使用无偏的自相关，会收获更加明显的周期性和整齐的峰值，但结果差别不大。

（三）
用延时相关估计距离是可行的，误差并不大。

d2的误差规律不明显，但是d1的测量值总是较d1的真实值要小，该方法有可能存在系统误差。

（四）
1.两信号频率相干图的峰值出现的位置也是两信号功率谱均具有峰值的频率。

2. 为了方便分析不同路信号相干图像，程序中选取128点为一段，并使用了主瓣更宽的blackman使图像更加平滑。

1）不同路同类型信号的相干图中可以观察到在功率谱中标注出δ(1-3Hz)、θ(4-7Hz)、α(8-13Hz)、β(14-30Hz)的范围均有对应峰值。

36Hz左右有未知原因的峰值， 46Hz后应该都是干扰性的噪声所以相干值保持很高。

这说明这些位置对应的频率分量大量包含在信号中。

2）可以看出与睁眼信号相比，闭眼信号除了δ波(1-3Hz)的值较低以外，θ(4-7Hz)、α(8-13Hz)、β(14-30Hz)和36Hz附近值均较高。

3.可以看出同路信号相干δ波(1-3Hz) ，β(14-30Hz)，36Hz附近有较明显峰值。

即睁闭眼信号频率成分的共性是δ波(1-3Hz) ，β(14-30Hz)，36Hz和噪声。

（五）
Matlab自带函数Xcorr()的定义和教科书不同，计算的是Ryx（m），并且计算出来的序列含有0，在Icorr中去掉了结果中由于x或有不等长产生的零。

十一、对本实验过程及方法、手段的改进建议：无
报告评分：
指导教师签字：。