Найдите сигнал в измерении

Вы получаете некоторые данные и хотели бы знать, совпадают ли они с более длинным потоком, вы измерились. Взаимная корреляция позволяет вам делать то определение, даже когда данные повреждаются шумом.

Загрузите в рабочую область запись звонка, вращающегося на поверхности стола. Обрежьте один второй фрагмент и слушайте его.

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

Смотрите также

Похожие темы