【生物医学信号处理及其MATLAB应用】Chapter 4 数字相关与数字卷积

4.1 线性相关

​ 线性相关是讨论两信号之间的同步性（Synchronism）或相似性（Similarity）或同相性（In-phase）或两信号的变化规律是否具有线性关系（Linear Relationship）或接近线性关系的程度。

​ 设有离散信号$x(n)$和$y(n)$，其线性相关函数为：

​ $r_{xy}(m)={\sum }_{n=-\infty }^{+\infty }x(n)y(n+m)$

​ 式中，$m$表示位移量，$m>0$表示$y(n)$序列左移，$m<0$表示右移，不同的$m$得到不同的$r_{xy}(m)$值。$r_{xy}(m)$值大于$0$表示有同相成分存在；小于$0$表示有反相成分存在；等于$0$表示两序列正交。线性相关运算的简洁表示为：

​ $r_{xy}(m)=x(n)\cdot y(n)$

​ 式中，“·”表示线性相关运算符（Correlation Operator）。

​ 当$x(n)$和$y(n)$完全相等时，就由互相关函数变成自相关函数。

4.2 循环相关

​ 设有序列$x(n)$和$y(n)$，其$N$点循环相关函数为：

​ $r_{xy}(m)={\sum }_{n=0}^{N-1}x(n)y((n+m){)}_N{R}_N(n)$

​ 循环相关运算的简洁表示为：

​ $r_{xy}(m)=x(n)⨀y(n)$

​ 当$x(n)$和$y(n)$完全相等时，由互相关函数编程自相关函数。

​ 在求循环相关时，序列要循环左移，一个$N$点长的序列$y(n)$行向量，进行$N$次循环移位得到一个移位矩阵，计算循环相关时，只要用另外一个行向量$x(n)$乘上该矩阵即可得到循环相关结果。下面给出$Matlab$中的循环左移函数$circlel()$:

function v = circlel(y)
% 输入向量必须是行向量
v = zeros(N, N);
for i = 1:N
	for j = 1:N
		v(i, j) = y(j);
	end
	L = y(1);
	for k = 1:N - 1
		y(k) = y(k + 1);
	end
	y(N) = L;
end

​ 给定行向量序列$x(n)$和$y(n)$，计算过程：$V=circlel(y);r=x\ast V;$即可。如果序列的长度不同，则把短的序列进行补零，使得两者点长相同，然后计算同上。

4.3 相干函数

​ 相关函数$r_{xy}(m)$值的大小与$x(n)$和$y(n)$序列值的大小有关，难以比较各组信号的相关程度，因而提出了相干函数（Coherent Function）的概念。

​ 时域相干函数：${\rho }_{xy}(m)=\frac{r_{x}y(m)}{[r_{xx}(0)\cdot r_{yy}(0){]}^{\frac{1}{2}}}$

​ 频域相干函数：${\gamma }_{xy}(k)=\frac{|P_{xy}(k){|}^2}{P_x(k)P_y(k)}$

4.4 线性卷积

​ 设有离散信号$x(n)$和$y(n)$，其线性卷积为：

​ $C_{xy}(m)={\sum }_{n=-\infty }^{+\infty }x(n)y(m-n)$

​ 与线性相关运算不同的是：

​ ①卷积运算时，$y(n)$要先反折得到$y(-n)$

​ ②$m>0$表示$y(-n)$序列右移，$m<0$表示左移，不同的$m$得到不同的$C_{xy}(m)$值。其余与相关运算相同。

​ 线性卷积运算的简洁表示为：

​ $C_{xy}(m)=x(n)\ast y(n)$

4.5 循环卷积

​ 有限长序列的循环移位是指$y((m-n){)}_N{R}_N(n)$，也就是先让序列$y(n)$以$N$为周期进行周期延拓，再进行反折，然后向右移位。只向一个方向移位的原因是：对周期序列向右移动一个位置，也就相当于向左移动了$N-1$个位置。最后取$(0,N-1)$的$N$个值就得到了循环移位后的$N$个序列值。

​ 设有序列$x(n)$和$y(n)$，其$N$点循环卷积为：

​ $C_{xy}(m)={\sum }_{n=0}^{N-1}x(n)y((m-n){)}_N{R}_N(n)$

4.6 MATLAB应用

4.6.1 成对相关

1. 在波形产生函数中选sinc波形发生器（s）和randn白噪声发生器（n），叠加成一个信号（x=s+weight*n），请改变白噪声信号的加权值，使得信噪比均为10，5，0，-5，-10dB，求得5个weight值，画出s和x信号，观察。

2. 用s信号与五个x信号，求他们的成对相关，corr或用corrcoef函数，画出相关系数与信噪比的关系图

3. 用s信号与x=weight*n-s做成对相关，画出相关系数与信噪比的关系图，观察s与x信号

clear;clc;close all
fs = 1000;  % 采样率设置
t = -10:1/fs:10;  % 设置时间变量
sn = sinc(t);  % 产生sinc信号
figure;subplot(3, 2, 1);
plot(t, sn);
title('sinc信号');xlabel('时间/s');  % 绘制sinc信号

wn = randn(1, length(sn));  % 产生与信号相同长度的随机白噪声
subplot(3, 2, 2);
plot(t, wn);
title('白噪信号');xlabel('时间/s');  % 绘制白噪信号

% 改变白噪声信号的加权值，设置信噪比为10，5，0，-5，-10dB，计算5个weight值
snr = [10 5 0 -5 -10];
w = sqrt(var(sn)./(var(wn).*(10.^(snr./10))));
% 获得五个观测信号x
x = [];
% 画出原始信号s和5个观测信号x的时域图，对比，观察
for i = 1:length(snr)
    x(:, i) = sn + w(i)*wn;
    str = ['信噪比', num2str(snr(i)), 'dB'];
    subplot(3, 2, i+1)
    plot(x(:, i));
    title(str);
    % 计算信号s和5个叠加信号x的成对相关系数
    corr_sx(i) = corr(sn', x(:, i));
end
figure;
plot(snr, corr_sx);
title('相关系数与信噪比的关系图');
xlabel('信噪比');ylabel('相关系数');

% 改变信噪比从-10dB到10dB，得到sinc信号与合成信号xn=sn-weight*wn的相关系数
snr1 = [-10 -5 0 5 10];
w1 = sqrt(var(sn)./(var(wn).*(10.^(snr1./10))));
figure;
x1 = [];
subplot(3, 2, 1);plot(t, sn);
title('sinc信号');xlabel('时间/s');
for i = 1:length(snr1)
    x1(:, i) = sn - w1(i)*wn;
    corr_sx1(i) = corr(sn', x1(:, i));
end

for i = 1:length(snr1)
    str = ['信噪比', num2str(snr1(i)), 'dB'];
    subplot(3, 2, i+1);
    plot(x1(:, i));
    title(str);
end

figure;
plot(snr1, corr_sx1);
title('相关系数与信噪比的关系图');
xlabel('信噪比');ylabel('相关系数');

运行结果：




Figure 1和Figure 3表明，信噪比越大，有用信号的特征在观测信号中表现得越不明显，即合成信号中的有用信号被噪声信号所淹没。Figure 2和Figure 4则对信噪比和相关系数的关系进行了计算，当$x=s+weight\ast n$时，原信号和加不同噪声的信号均为正相关，信噪比越大，相关系数也越大，而当$x=-s+weight\ast n$时，相关系数为负，信噪比越大，相关系数的绝对值越大。

4.6.2 相关——信号的检测

1. 任意产生一个白噪信号（w）和一个双频率的信号（s）

2. 信号s的周期是多少？仿真s，w和x，（x=s+a*w，a=1），计算三个信号的线性相关，并画出时间信号和线性相关结果

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

clear;clc;close all
fs = 100;  % 采样率1000Hz
n = -10:1/fs:10;
wn = randn(1, length(n));  % 构造白噪信号wn
figure;
subplot(3, 1, 1);
plot(n, wn);
title('白噪信号wn');
xlabel('时间/s');ylabel('幅度');
% 构造双频率的信号sn，计算出该信号的周期
sn = 2*sin(2*pi*n) + 4*sin(3*pi*n); 
subplot(3, 1, 2);
plot(n, sn);
title('双频率信号sn');
xlabel('时间/s');ylabel('幅度');
% 构造含噪双频信号xn（x=s+a*w，a=1）
xn = sn + wn;
subplot(3, 1, 3);
plot(n, xn);
title('含噪双频信号xn');
xlabel('时间/s');ylabel('幅度');

% 对双频信号求自相关
[corr_sn, lag_sn] = xcorr(sn, 'unbiased');
% 对含噪双频信号求自相关
[corr_xn, lag_xn] = xcorr(xn, 'unbiased');
figure;
subplot(2, 1, 1);
plot(lag_sn/fs, corr_sn);
title('双频信号的自相关函数');
subplot(2, 1, 2);
plot(lag_xn/fs, corr_xn);
title('含噪双频信号的自相关函数');

% 使用findpeaks函数估计含噪双频信号的周期大小
[pks, locs] = findpeaks(corr_xn, 'MinPeakDistance', 0.8, 'MinPeakHeight', 9);
figure;
plot(corr_xn);hold on
plot(locs, pks, '*');
T = diff(locs/fs);
T_mean = mean(T);

% 改变信噪比从-20dB到20dB，得到双频信号与合成信号x=s+a*w的信噪比-相关系数图
snr = [-20 -10 0 10 20];
w = sqrt(var(sn)./(var(wn).*(10.^(snr./10))));
% 获得五个观测信号x
x = [];
% 画出原始信号s和5个观测信号x的时域图，对比，观察
figure;subplot(3, 2, 1);
plot(n, sn);
title('原始信号');xlabel('时间/s');
for i = 1:length(snr)
    x(:, i) = sn + w(i)*wn;
    str = ['信噪比', num2str(snr(i)), 'dB'];
    subplot(3, 2, i+1);
    plot(x(:, i));
    title(str);
    % 计算信号s和5个叠加信号x的成对相关系数
    corr_sx(i) = corr(sn', x(:, i));
end
figure;
plot(snr, corr_sx);
title('相关系数与信噪比的关系图');
xlabel('信噪比');ylabel('相关系数');

for j = 1:length(snr)
    x1(:, j) = sn + w(j)*wn;
    str = ['信噪比', num2str(snr(j)), 'dB'];
    subplot(3, 2, 1);
    plot(n, sn);
    title('原始信号');xlabel('时间/s');
    subplot(3, 2, j+1);
    plot(xcorr(x1(:, j), 'unbiased'));
    title(str);
end

period1 = signal_period(x(:, 1), fs, 150, 50);
period2 = signal_period(x(:, 2), fs, 150, 9);
period3 = signal_period(x(:, 3), fs, 150, 9);
period4 = signal_period(x(:, 4), fs, 150, 9);
period5 = signal_period(x(:, 5), fs, 150, 9);
period = [period1 period2 period3 period4 period5];

figure;
plot(snr, period, 'o-');
xlabel('信噪比');ylabel('估计的周期');
title('信噪比-估计的周期');

以上代码调用了下述自编函数：

function T = signal_period(signal, Fs, min_distance, min_peak)
% signal: 输入信号向量
% Fs: 采样率
% T: 信号周期

% 计算自相关函数
R = xcorr(signal);

% 找到第一个峰值
[~, locs] = findpeaks(R, 'MinPeakDistance', min_distance, 'MinPeakHeight', min_peak);

% 计算周期
T = mean(diff(locs)/Fs);
end

运行结果：







采用双频信号$sn=2sin(2\pi n)+4sin(3\pi n)$，计算可知双频信号的周期为$2$。在对信噪比-估计周期值关系进行分析前，先使用$findpeaks$函数对含噪双频信号的周期大小进行了估计，此处取$MinPeakDistance=1.9$，$MinPeakHeight=9$，计算得到周期值为$1.7555$，与实际观测结果间误差较小，基本相符。该条结论说明我们可以从线性相关函数中检测出有用信号的周期，即因为 ，因此得到自相关运算的计算结果后，对比自相关序列与原始周期序列的周期是否相同，即可寻找被噪声淹没的周期信号的周期值。

4.6.3 延时相关——估计距离的一种方法

clear;clc;close 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，d2两个距离的参数
d1 = 340;
d2 = 680;
t1 = d1/c0;
t2 = (d1 + 2*d2)/c0;
tc = 2*(d1 + d2)/c0;
Lag = 500;  % 最大延迟样点数
pt = sinc(2*pi*t);  % 原信号
figure;plot(t, pt);
xlabel('时间/s');title('原始信号p(t)');

% 画出一号传声器输出x(t)
xt = pt + A1*sinc(2*pi*(t - tc));
figure;subplot(2, 1, 1);
plot(t, xt);
title('一号传声器输出x(t)信号');
% 画出二号传声器输出y(t)
yt = A*sinc(2*pi*(t - t1)) + A2*sinc(2*pi*(t - t2));
subplot(2, 1, 2);
plot(t, yt);
title('二号传声器输出y(t)信号');

% 求三个信号的自相关及x与y信号的互相关
[Rpp, lags] = xcorr(pt, 'biased');
[Rxx, lagx] = xcorr(xt, 'biased');
[Ryy, lagt] = xcorr(yt, 'biased');
[Rxy, lagy] = xcorr(xt, yt, 'biased');
% 画出以上四个相关函数图
figure;
subplot(2, 2, 1);rts = lags/Fs;
plot(rts, Rpp);
title('p(t)信号自相关');
subplot(2, 2, 2);rts = lagx/Fs;
plot(rts, Rxx);
title('x(t)信号自相关');
subplot(2, 2, 3);rts = lagt/Fs;
plot(rts, Ryy);
title('y(t)信号自相关');
subplot(2, 2, 4);rts = lagy/Fs;
plot(rts, Rxy);
title('x(t)、y(t)信号互相关');

% 找互相关函数的四个极大值点对应横坐标时间（findpeaks）
r_wave1 = [];
Threshold = (max(Rxy) - min(Rxy))*0.5 + min(Rxy);
[r_wave1(:,1), r_wave1(:,2)] = findpeaks(Rxy, 50, 'MinPeakDistance', 0.025, 'MinPeakHeight', Threshold);
r_wave1(:,2) = r_wave1(:, 2) + min(rts);
figure;
plot(rts, Rxy);
hold on; 
plot(r_wave1(:,2), r_wave1(:,1), '*');
title('x(t)、y(t)信号互相关和距离检测');
% 计算d1和d2
get_d1 = c0*r_wave1(3, 2)
get_d2 = 0.5*(c0*r_wave1(4, 2) - get_d1)