Вы получаете некоторые данные и хотели бы знать, соответствуют ли они более длинному потоку данных, которые вы измерили. Взаимная корреляция позволяет вам сделать такое определение, даже если данные искажены шумом.
Загрузите в рабочее пространство запись вращения кольца на столе. Обрежьте односекундный фрагмент и прослушайте его.
load('Ring.mat')
Time = 0:1/Fs:(length(y)-1)/Fs;
m = min(y);
M = max(y);
Full_sig = double(y);
timeA = 7;
timeB = 8;
snip = timeA*Fs:timeB*Fs;
Fragment = Full_sig(snip);
% To hear, type soundsc(Fragment,Fs)
Постройте график сигнала и фрагмента. Выделите конечные точки фрагмента для референса.
plot(Time,Full_sig,[timeA timeB;timeA timeB],[m m;M M],'r--')
xlabel('Time (s)')
ylabel('Clean')
axis tight
plot(snip/Fs,Fragment)
xlabel('Time (s)')
ylabel('Clean')
title('Fragment')
axis tight
Вычислите и постройте график взаимной корреляции полного сигнала и фрагмента.
[xCorr,lags] = xcorr(Full_sig,Fragment);
plot(lags/Fs,xCorr)
grid
xlabel('Lags (s)')
ylabel('Clean')
axis tight
Задержка, при которой взаимная корреляция является наибольшей, - это временная задержка между начальными точками сигналов. Переместите сигнал и наложите фрагмент.
[~,I] = max(abs(xCorr));
maxt = lags(I);
Trial = NaN(size(Full_sig));
Trial(maxt+1:maxt+length(Fragment)) = Fragment;
plot(Time,Full_sig,Time,Trial)
xlabel('Time (s)')
ylabel('Clean')
axis tight
Повторите процедуру, но добавьте шум отдельно к сигналу и фрагменту. Звук невозможно выделить из шума.
NoiseAmp = 0.2*max(abs(Fragment));
Fragment = Fragment+NoiseAmp*randn(size(Fragment));
Full_sig = Full_sig+NoiseAmp*randn(size(Full_sig));
% To hear, type soundsc(Fragment,Fs)
plot(Time,Full_sig,[timeA timeB;timeA timeB],[m m;M M],'r--')
xlabel('Time (s)')
ylabel('Noisy')
axis tight
Процедура находит недостающий фрагмент, несмотря на высокий уровень шума.
[xCorr,lags] = xcorr(Full_sig,Fragment);
plot(lags/Fs,xCorr)
grid
xlabel('Lags (s)')
ylabel('Noisy')
axis tight
[~,I] = max(abs(xCorr));
maxt = lags(I);
Trial = NaN(size(Full_sig));
Trial(maxt+1:maxt+length(Fragment)) = Fragment;
figure
plot(Time,Full_sig,Time,Trial)
xlabel('Time (s)')
ylabel('Noisy')
axis tight