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

电子科技大学生命科学与技术学院标准实验报告

（实验）课程名称生物医学信号处理

2018-2019-第2学期

电子科技大学教务处制表

一、实验室名称：品学楼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.线性相关

设有离散信号 、 ，则其线性相关函数为

[]0

0(0)(1)(1)0(0)

(1)(2)

0(0)(1)

(1)(1)

(2)

00(1)

(0)

(2)(1)0

0y 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 结果进行比较。

1

()()()

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互相关图

表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)