Matlab_xcorr_互相关函数的讨论

Matlab_xcorr_互相关函数的讨论
假设两个平稳信号x和y，如果x(t+τ)=y(t) ，则可通过互相关求τ。

由于信号处理相关知识都蘸酱吃了，理论相关的部分咱们来⽇⽅长（我⼀定可能会补充的）。

XCORR 实现
⾸先，通过实现 xcorr 函数介绍：
clc
clear
close
% 实现 xcorr 函数
% 基本设置
T = 1; % [s] 总时间长度
fs = 5000; % [Hz] 采样频率
t = 0:1/fs:T; % [s] 时间坐标
N = length(t); % 信号个数
% 信号⽣成
tm = [ t(1:N) - T , t(2:N) ]; % 相关结果的时间延迟坐标轴
td1 = 0.2*T; % x 信号时间延迟
td2 = 0.3*T; % y 信号时间延迟
noise = rand(1,2*N); % ⽣成了两倍时间 T 长度的噪声 [0,1]噪声
x = noise(1+round(td1*fs):N+round(td1*fs))-0.5*ones(1,N);
y = noise(1+round(td2*fs):N+round(td2*fs))-0.5*ones(1,N);
% 求取互相关
z1 = xcorr(x,y); % Matlab ⾃带函数
[~,I1] = max(abs(z1));　% 模仿 Matlab doc 给出延迟坐标
z2 = zeros(1,N); % ⾃编函数
for n = 1:length(tm)
z2(n) = sum( x( max(1,n-N+1):min(n,N) ).*y( max(1,N-n+1):min(2*N-n,N) ) );
end
[~,I2] = max(abs(z2));　% 模仿 Matlab doc 给出延迟坐标
%---------------------计算说明--------------------%
% case1： | case2: %
% .N | .2*N-n %
% y: .......... | y: .......... %
% .N-n+1 | .1 %
% .n | .N %
% x: .......... | x: .......... %
% .1 | .n-N+1 %
%------------------------------------------------%
err = z1-z2; % 两种算法的差
% 绘图
subplot(1,3,1)
plot(tm,z1)
title('Matlab function')
xlabel('time delay')
ylabel('Amp')
a1 = gca;
a1.XTick = sort([-1:0.5:1 tm(I1)]);
subplot(1,3,2)
plot(tm,z2)
title('My function')
xlabel('time delay')
ylabel('Amp')
a2 = gca;
a2.XTick = sort([-1:0.5:1 tm(I2)]);
subplot(1,3,3)
plot(tm,err,'.-')
title('error')
xlabel('time delay')
ylabel('Amp')
suptitle('xcorr realization')
以上 Matlab 代码可以得到下⾯的结果。

从左到右依次是 Matlab ⾃带函数、我编的互相关函数、两个函数的差值。

不难发现：两个函数⼗分接近，但是差值不为零。

个⼈猜测是因为 xcorr 的求和和 sum 求和的截断误差不同所致。

这个误差的来源我懒得去编程序验证了——毕竟10-16量级的差别，没多⼤深究的意义。

但是可以注意到这个差值有四个特点：
- ⼩幅值时有固定⼏个数值
- 每跑⼀次程序，rand 产⽣的噪声数据不同，error 值不同
- 呈“纺锤型”，中间⾼，两边低
- 实际值⼤的数据点，error 值⼤
最后要谈⼀下 xcorr 的噪声问题。

我们通常使⽤的噪声是⽩噪声，或者⾼斯⽩，有⼀个很重要的特点就是均值为零，也就是说没有直流分量。

但是当我们的噪声存在直流分量的时候（⽐如上⾯的噪声信号直接使⽤rand(1,2*N)时），互相关就是⼀个类似等腰三⾓形的东西了（想想门函数卷积）。

回忆⼀下，对于存在稳定周期分量的两组信号x、y⽽⾔，互相关结果将会是⼀个幅度为“纺锤形”的周期震荡的信号。

由此可观：互相关⼀⽅⾯可以得到⾮周期信号延迟结果，同时也能反映极端情况下，相同频率成分的存在，这⼀点可以⽤来观察⼯频⼲扰程度。

XCORR 与 CONV
互相关 xcorr 与 conv 的差别在于两点：
- xcorr 在两段信号较短者后补零，使两段信号长度⼀致
- xcorr 直接⽤两个信号的各种延迟做相乘求和，conv 使⽤翻褶后的信号做相乘求和
这导致了：
1、xcorr(x,y) 中 (x,y) 顺序有影响，⽽conv(x,y) 没有
2、两者在⼤部分情况下得到的结果是不⼀样的，但是对于⼀些有趣的对称信号是存在等价关系的。

有兴趣的读者可以搞⼀搞，找找规律。

因为本⼈并不搞对称相关的研究，这点就不展开了。

下⾯的例⼦是有等价关系的。

clc
clear
close
% ⽐较 conv xcorr
% 例⼦
A = ones(1,12); % -3:3
B = 0:4; % 3:-1:-3
C = xcorr(A,B);
D = conv(A,B);
%绘图
subplot(2,2,1)
plot(A,'.-')
ylim([ -0.1 5.1 ])
xlim([ 0.9 12.1])
title('A = ones(1,12)')
xlabel('n')
ylabel('Amp')
subplot(2,2,2)
plot(B,'.-')
ylim([ -0.1 5.1 ])
xlim([ 0.9 12.1])
title('B = 0:4')
xlabel('n')
ylabel('Amp')
subplot(2,2,3)
plot(C,'.-')
ylim([ -0.1 15.1 ])
xlim([ 0.9 25.1])
title('xcorr 结果')
xlabel('n')
ylabel('Amp')
subplot(2,2,4)
plot(D,'.-')
ylim([ -0.1 15.1 ])
xlim([ 0.9 25.1])
title('cone 结果')
xlabel('n')
ylabel('Amp')
suptitle('conv与xcorr对⽐')
有兴趣的读者可以试着⽤给定函数实现⽬标函数： - xcorr --> fliplr
- xcorr --> conv
- conv --> fliplr
- conv --> xcorr
END
Processing math: 100%。