Найдите сигнал в измерении
Вы получаете некоторые данные и хотели бы знать, соответствует ли они более длинному потоку, который вы измерили. Перекрестная корреляция позволяет вам сделать это определение, даже когда данные повреждены шумом.
Загрузите в рабочую область запись звонка, вращающейся на столешнице. Обрезать фрагмент на одну секунду и послушать его.
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
