matlab的xcorr函数,Matlab_xcorr_互相关函数的讨论
假設(shè)兩個(gè)平穩(wěn)信號(hào) $\textbf{x}$ 和 $\textbf{y}$ ,如果 $x\left(t+\tau\right)= y\left(t\right)$ ,則可通過(guò)互相關(guān)求 $\tau$ 。由于信號(hào)處理相關(guān)知識(shí)都蘸醬吃了,理論相關(guān)的部分咱們來(lái)日方長(zhǎng)(我一定可能會(huì)補(bǔ)充的)。
XCORR 實(shí)現(xiàn)
首先,通過(guò)實(shí)現(xiàn) xcorr 函數(shù)介紹互相關(guān)計(jì)算流程:
clc
clear
close
% 實(shí)現(xiàn) xcorr 函數(shù)
% 基本設(shè)置
T = 1; % [s] 總時(shí)間長(zhǎng)度
fs = 5000; % [Hz] 采樣頻率
t = 0:1/fs:T; % [s] 時(shí)間坐標(biāo)
N = length(t); % 信號(hào)個(gè)數(shù)
% 信號(hào)生成
tm = [ t(1:N) - T , t(2:N) ]; % 相關(guān)結(jié)果的時(shí)間延遲坐標(biāo)軸
td1 = 0.2*T; % x 信號(hào)時(shí)間延遲
td2 = 0.3*T; % y 信號(hào)時(shí)間延遲
noise = rand(1,2*N); % 生成了兩倍時(shí)間 T 長(zhǎng)度的噪聲 [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);
% 求取互相關(guān)
z1 = xcorr(x,y); % Matlab 自帶函數(shù)
[~,I1] = max(abs(z1)); % 模仿 Matlab doc 給出延遲坐標(biāo)
z2 = zeros(1,N); % 自編函數(shù)
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 給出延遲坐標(biāo)
%---------------------計(jì)算說(shuō)明--------------------%
% 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 代碼可以得到下面的結(jié)果。從左到右依次是 Matlab 自帶函數(shù)、我編的互相關(guān)函數(shù)、兩個(gè)函數(shù)的差值。不難發(fā)現(xiàn):兩個(gè)函數(shù)十分接近,但是差值不為零。個(gè)人猜測(cè)是因?yàn)?xcorr 的求和和 sum 求和的截?cái)嗾`差不同所致。這個(gè)誤差的來(lái)源我懶得去編程序驗(yàn)證了——畢竟10-16量級(jí)的差別,沒(méi)多大深究的意義。但是可以注意到這個(gè)差值有四個(gè)特點(diǎn):
- 小幅值時(shí)有固定幾個(gè)數(shù)值
- 每跑一次程序,rand 產(chǎn)生的噪聲數(shù)據(jù)不同,error 值不同
- 呈“紡錘型”,中間高,兩邊低
- 實(shí)際值大的數(shù)據(jù)點(diǎn),error 值大
最后要談一下 xcorr 的噪聲問(wèn)題。我們通常使用的噪聲是白噪聲,或者高斯白,有一個(gè)很重要的特點(diǎn)就是均值為零,也就是說(shuō)沒(méi)有直流分量。但是當(dāng)我們的噪聲存在直流分量的時(shí)候(比如上面的噪聲信號(hào)直接使用rand(1,2*N)時(shí)),互相關(guān)就是一個(gè)類似等腰三角形的東西了(想想門函數(shù)卷積)。回憶一下,對(duì)于存在穩(wěn)定周期分量的兩組信號(hào) $\textbf{x}$ 、 $\textbf{y}$ 而言,互相關(guān)結(jié)果將會(huì)是一個(gè)幅度為“紡錘形”的周期震蕩的信號(hào)。由此可觀:互相關(guān)一方面可以得到非周期信號(hào)延遲結(jié)果,同時(shí)也能反映極端情況下,相同頻率成分的存在,這一點(diǎn)可以用來(lái)觀察工頻干擾程度。
XCORR 與 CONV
互相關(guān) xcorr 與 conv 的差別在于兩點(diǎn):
- xcorr 在兩段信號(hào)較短者后補(bǔ)零,使兩段信號(hào)長(zhǎng)度一致
- xcorr 直接用兩個(gè)信號(hào)的各種延遲做相乘求和,conv 使用翻褶后的信號(hào)做相乘求和
這導(dǎo)致了:
1、xcorr(x,y) 中 (x,y) 順序有影響,而conv(x,y) 沒(méi)有
2、兩者在大部分情況下得到的結(jié)果是不一樣的,但是對(duì)于一些有趣的對(duì)稱信號(hào)是存在等價(jià)關(guān)系的。有興趣的讀者可以搞一搞,找找規(guī)律。因?yàn)楸救瞬⒉桓銓?duì)稱相關(guān)的研究,這點(diǎn)就不展開(kāi)了。下面的例子是有等價(jià)關(guān)系的。
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 結(jié)果')
xlabel('n')
ylabel('Amp')
subplot(2,2,4)
plot(D,'.-')
ylim([ -0.1 15.1 ])
xlim([ 0.9 25.1])
title('cone 結(jié)果')
xlabel('n')
ylabel('Amp')
suptitle('conv與xcorr對(duì)比')
有興趣的讀者可以試著用給定函數(shù)實(shí)現(xiàn)目標(biāo)函數(shù):
- xcorr --> fliplr
- xcorr --> conv
- conv --> fliplr
- conv --> xcorr
END
總結(jié)
以上是生活随笔為你收集整理的matlab的xcorr函数,Matlab_xcorr_互相关函数的讨论的全部?jī)?nèi)容,希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇: 安装反射大师具体步骤与使用教程
- 下一篇: Android OKHTTP发起请求提示