Chapter 4 数字相关与数字卷积
4.1 线性相关
线性相关是讨论两信号之间的同步性(Synchronism)或相似性(Similarity)或同相性(In-phase)或两信号的变化规律是否具有线性关系(Linear Relationship)或接近线性关系的程度。
设有离散信号 x ( n ) x(n) x(n)和 y ( n ) y(n) y(n),其线性相关函数为:
r x y ( m ) = ∑ n = − ∞ + ∞ x ( n ) y ( n + m ) r_{xy}(m)=\sum_{n=-∞}^{+∞}x(n)y(n+m) rxy(m)=∑n=−∞+∞x(n)y(n+m)
式中, m m m表示位移量, m > 0 m>0 m>0表示 y ( n ) y(n) y(n)序列左移, m < 0 m<0 m<0表示右移,不同的 m m m得到不同的 r x y ( m ) r_{xy}(m) rxy(m)值。 r x y ( m ) r_{xy}(m) rxy(m)值大于 0 0 0表示有同相成分存在;小于 0 0 0表示有反相成分存在;等于 0 0 0表示两序列正交。线性相关运算的简洁表示为:
r x y ( m ) = x ( n ) ⋅ y ( n ) r_{xy}(m)=x(n)·y(n) rxy(m)=x(n)⋅y(n)
式中,“·”表示线性相关运算符(Correlation Operator)。
当 x ( n ) x(n) x(n)和 y ( n ) y(n) y(n)完全相等时,就由互相关函数变成自相关函数。
4.2 循环相关
设有序列 x ( n ) x(n) x(n)和 y ( n ) y(n) y(n),其 N N N点循环相关函数为:
r x y ( m ) = ∑ n = 0 N − 1 x ( n ) y ( ( n + m ) ) N R N ( n ) r_{xy}(m)=\sum_{n=0}^{N-1}x(n)y((n+m))_{N}R_N(n) rxy(m)=∑n=0N−1x(n)y((n+m))NRN(n)
循环相关运算的简洁表示为:
r x y ( m ) = x ( n ) ⨀ y ( n ) r_{xy}(m)=x(n)\bigodot y(n) rxy(m)=x(n)⨀y(n)
当 x ( n ) x(n) x(n)和 y ( n ) y(n) y(n)完全相等时,由互相关函数编程自相关函数。
在求循环相关时,序列要循环左移,一个 N N N点长的序列 y ( n ) y(n) y(n)行向量,进行 N N N次循环移位得到一个移位矩阵,计算循环相关时,只要用另外一个行向量 x ( n ) x(n) x(n)乘上该矩阵即可得到循环相关结果。下面给出 M a t l a b Matlab Matlab中的循环左移函数 c i r c l e l ( ) circlel( ) 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 ) x(n) x(n)和 y ( n ) y(n) y(n),计算过程: V = c i r c l e l ( y ) ; r = x ∗ V ; V=circlel(y);r=x*V; V=circlel(y);r=x∗V;即可。如果序列的长度不同,则把短的序列进行补零,使得两者点长相同,然后计算同上。
4.3 相干函数
相关函数 r x y ( m ) r_{xy}(m) rxy(m)值的大小与 x ( n ) x(n) x(n)和 y ( n ) y(n) y(n)序列值的大小有关,难以比较各组信号的相关程度,因而提出了相干函数(Coherent Function)的概念。
时域相干函数: ρ x y ( m ) = r x y ( m ) [ r x x ( 0 ) ⋅ r y y ( 0 ) ] 1 2 ρ_{xy}(m)=\frac{r_xy(m)}{[r_{xx}(0)·r_{yy}(0)]^{\frac{1}{2}}} ρxy(m)=[rxx(0)⋅ryy(0)]21rxy(m)
频域相干函数: γ x y ( k ) = ∣ P x y ( k ) ∣ 2 P x ( k ) P y ( k ) \gamma_{xy}(k)=\frac{|P_{xy}(k)|^2}{P_x(k)P_y(k)} γxy(k)=Px(k)Py(k)∣Pxy(k)∣2
4.4 线性卷积
设有离散信号 x ( n ) x(n) x(n)和 y ( n ) y(n) y(n),其线性卷积为:
C x y ( m ) = ∑ n = − ∞ + ∞ x ( n ) y ( m − n ) C_{xy}(m)=\sum_{n=-∞}^{+∞}x(n)y(m-n) Cxy(m)=∑n=−∞+∞x(n)y(m−n)
与线性相关运算不同的是:
①卷积运算时, y ( n ) y(n) y(n)要先反折得到 y ( − n ) y(-n) y(−n)
② m > 0 m>0 m>0表示 y ( − n ) y(-n) y(−n)序列右移, m < 0 m<0 m<0表示左移,不同的 m m m得到不同的 C x y ( m ) C_{xy}(m) Cxy(m)值。其余与相关运算相同。
线性卷积运算的简洁表示为:
C x y ( m ) = x ( n ) ∗ y ( n ) C_{xy}(m)=x(n)*y(n) Cxy(m)=x(n)∗y(n)
4.5 循环卷积
有限长序列的循环移位是指 y ( ( m − n ) ) N R N ( n ) y((m-n))_{N}R_N(n) y((m−n))NRN(n),也就是先让序列 y ( n ) y(n) y(n)以 N N N为周期进行周期延拓,再进行反折,然后向右移位。只向一个方向移位的原因是:对周期序列向右移动一个位置,也就相当于向左移动了 N − 1 N-1 N−1个位置。最后取 ( 0 , N − 1 ) (0,N-1) (0,N−1)的 N N N个值就得到了循环移位后的 N N N个序列值。
设有序列 x ( n ) x(n) x(n)和 y ( n ) y(n) y(n),其 N N N点循环卷积为:
C x y ( m ) = ∑ n = 0 N − 1 x ( n ) y ( ( m − n ) ) N R N ( n ) C_{xy}(m)=\sum_{n=0}^{N-1}x(n)y((m-n))_{N}R_N(n) Cxy(m)=∑n=0N−1x(n)y((m−n))NRN(n)
4.6 MATLAB应用
4.6.1 成对相关
-
在波形产生函数中选sinc波形发生器(s)和randn白噪声发生器(n),叠加成一个信号(x=s+weight*n),请改变白噪声信号的加权值,使得信噪比均为10,5,0,-5,-10dB,求得5个weight值,画出s和x信号,观察。
-
用s信号与五个x信号,求他们的成对相关,corr或用corrcoef函数,画出相关系数与信噪比的关系图
-
用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
+
w
e
i
g
h
t
∗
n
x=s+weight*n
x=s+weight∗n时,原信号和加不同噪声的信号均为正相关,信噪比越大,相关系数也越大,而当
x
=
−
s
+
w
e
i
g
h
t
∗
n
x=-s+weight*n
x=−s+weight∗n时,相关系数为负,信噪比越大,相关系数的绝对值越大。
4.6.2 相关——信号的检测
-
任意产生一个白噪信号(w)和一个双频率的信号(s)
-
信号s的周期是多少?仿真s,w和x,(x=s+a*w,a=1),计算三个信号的线性相关,并画出时间信号和线性相关结果
-
改变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
运行结果:
采用双频信号
s
n
=
2
s
i
n
(
2
π
n
)
+
4
s
i
n
(
3
π
n
)
sn=2sin(2\pi n)+4sin(3\pi n)
sn=2sin(2πn)+4sin(3πn),计算可知双频信号的周期为
2
2
2。在对信噪比-估计周期值关系进行分析前,先使用
f
i
n
d
p
e
a
k
s
findpeaks
findpeaks函数对含噪双频信号的周期大小进行了估计,此处取
M
i
n
P
e
a
k
D
i
s
t
a
n
c
e
=
1.9
MinPeakDistance=1.9
MinPeakDistance=1.9,
M
i
n
P
e
a
k
H
e
i
g
h
t
=
9
MinPeakHeight=9
MinPeakHeight=9,计算得到周期值为
1.7555
1.7555
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)