Сравнение сигналов с различными скоростями дискретизации

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

load relatedsig.mat

figure
ax(1) = subplot(3,1,1);
plot((0:numel(T1)-1)/Fs1,T1,'k')
ylabel('Template 1')
grid on
ax(2) = subplot(3,1,2); 
plot((0:numel(T2)-1)/Fs2,T2,'r')
ylabel('Template 2')
grid on
ax(3) = subplot(3,1,3); 
plot((0:numel(S)-1)/Fs,S)
ylabel('Signal')
grid on
xlabel('Time (secs)')
linkaxes(ax(1:3),'x')
axis([0 1.61 -4 4])

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

[Fs1 Fs2 Fs]

ans = 1×3

        4096        4096        8192

Различные длины мешают вычислить разницу между двумя сигналами, но это можно легко исправить, извлекая общую часть сигналов. Кроме того, не всегда необходимо выравнивать длины. Перекрестная корреляция может выполняться между сигналами различной длины, но важно обеспечить, чтобы они имели идентичные частоты дискретизации. Самым безопасным способом для этого является повторная выборка сигнала с более низкой частотой дискретизации. resample функция применяет сглаживающий (низкочастотный) фильтр FIR к сигналу во время процесса повторной дискретизации.

[P1,Q1] = rat(Fs/Fs1);          % Rational fraction approximation
[P2,Q2] = rat(Fs/Fs2);          % Rational fraction approximation
T1 = resample(T1,P1,Q1);        % Change sampling rate by rational factor
T2 = resample(T2,P2,Q2);        % Change sampling rate by rational factor

Поиск сигнала в измерении

Теперь мы можем перекрестно коррелировать сигнал S с шаблонами T1 и T2 с xcorr для определения наличия совпадения.

[C1,lag1] = xcorr(T1,S);        
[C2,lag2] = xcorr(T2,S);        

figure
ax(1) = subplot(2,1,1); 
plot(lag1/Fs,C1,'k')
ylabel('Amplitude')
grid on
title('Cross-correlation between Template 1 and Signal')
ax(2) = subplot(2,1,2); 
plot(lag2/Fs,C2,'r')
ylabel('Amplitude') 
grid on
title('Cross-correlation between Template 2 and Signal')
xlabel('Time(secs)') 
axis(ax(1:2),[-1.5 1.5 -700 700 ])

Первый субплот указывает, что сигнал и шаблон 1 менее коррелированы, в то время как высокий пик во втором субплоте указывает, что сигнал присутствует во втором шаблоне.

[~,I] = max(abs(C2));
SampleDiff = lag2(I)

SampleDiff = 499

timeDiff = SampleDiff/Fs

timeDiff = 0.0609

Пик взаимной корреляции подразумевает, что сигнал присутствует в шаблоне T2 начиная с 61 мс. Другими словами, сигнал T2 ведет сигнал S на 499 выборок, как указано SampleDiff. Эта информация может использоваться для выравнивания сигналов.

Измерение задержки между сигналами и их выравнивание

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

figure
ax(1) = subplot(3,1,1);
plot(s1)
ylabel('s1')
grid on
ax(2) = subplot(3,1,2); 
plot(s2,'k')
ylabel('s2')
grid on
ax(3) = subplot(3,1,3); 
plot(s3,'r')
ylabel('s3')
grid on
xlabel('Samples')
linkaxes(ax,'xy')

Мы также можем использовать finddelay функция для поиска задержки между двумя сигналами.

t21 = finddelay(s1,s2)

t21 = -350

t31 = finddelay(s1,s3)

t31 = 150

t21 указывает, что s2 отстает s1 на 350 выборок, а t31 указывает, что s3 ведет s1 на 150 выборок. Эта информация теперь может использоваться для выравнивания 3 сигналов путем временного сдвига сигналов. Мы также можем использовать alignsignals непосредственно выравнивать сигналы, которые выравнивают два сигнала, задерживая самый ранний сигнал.

s1 = alignsignals(s1,s3);
s2 = alignsignals(s2,s3);

figure
ax(1) = subplot(3,1,1);
plot(s1)
grid on 
title('s1')
axis tight
ax(2) = subplot(3,1,2);
plot(s2)
grid on 
title('s2')
axis tight
ax(3) = subplot(3,1,3); 
plot(s3)
grid on 
title('s3')
axis tight
linkaxes(ax,'xy')

Сравнение частотного содержания сигналов

Спектр мощности отображает мощность, присутствующую на каждой частоте. Спектральная когерентность идентифицирует корреляцию частотной области между сигналами. Значения когерентности, имеющие тенденцию к 0, указывают, что соответствующие частотные компоненты не коррелированы, в то время как значения, имеющие тенденцию к 1, указывают, что соответствующие частотные компоненты коррелированы. Рассмотрим два сигнала и их соответствующие спектры мощности.

Fs = FsSig;         % Sampling Rate

[P1,f1] = periodogram(sig1,[],[],Fs,'power');
[P2,f2] = periodogram(sig2,[],[],Fs,'power');

figure
t = (0:numel(sig1)-1)/Fs;
subplot(2,2,1)
plot(t,sig1,'k')
ylabel('s1')
grid on
title('Time Series')
subplot(2,2,3)
plot(t,sig2)
ylabel('s2')
grid on
xlabel('Time (secs)')
subplot(2,2,2)
plot(f1,P1,'k')
ylabel('P1')
grid on
axis tight
title('Power Spectrum')
subplot(2,2,4)
plot(f2,P2)
ylabel('P2')
grid on
axis tight
xlabel('Frequency (Hz)')

mscohere функция вычисляет спектральную когерентность между двумя сигналами. Это подтверждает, что sig1 и sig2 имеют два коррелированных компонента около 35 Гц и 165 Гц. В частотах, где спектральная когерентность высока, относительная фаза между коррелированными компонентами может быть оценена с фазой перекрестного спектра.

[Cxy,f] = mscohere(sig1,sig2,[],[],[],Fs);
Pxy = cpsd(sig1,sig2,[],[],[],Fs);
phase = -angle(Pxy)/pi*180;
[pks,locs] = findpeaks(Cxy,'MinPeakHeight',0.75);

figure
subplot(2,1,1)
plot(f,Cxy)
title('Coherence Estimate')
grid on
hgca = gca;
hgca.XTick = f(locs);
hgca.YTick = 0.75;
axis([0 200 0 1])
subplot(2,1,2)
plot(f,phase)
title('Cross-spectrum Phase (deg)')
grid on
hgca = gca;
hgca.XTick = f(locs); 
hgca.YTick = round(phase(locs));
xlabel('Frequency (Hz)')
axis([0 200 -180 180])

Фазовое отставание между компонентами 35 Гц близко к -90 градусам, а фазовое отставание между компонентами 165 Гц близко к -60 градусам.

Поиск периодичности в сигнале

Рассмотрим набор измерений температуры в офисном здании в зимний сезон. Измерения проводили каждые 30 минут в течение примерно 16,5 недель.

load officetemp.mat  

Fs = 1/(60*30);                 % Sample rate is 1 sample every 30 minutes
days = (0:length(temp)-1)/(Fs*60*60*24); 

figure
plot(days,temp)
title('Temperature Data')
xlabel('Time (days)')
ylabel('Temperature (Fahrenheit)')
grid on

С температурами в низких 70-х нужно убрать среднее для анализа небольших колебаний сигнала. xcov функция удаляет среднее значение сигнала перед вычислением взаимной корреляции. Он возвращает перекрестную ковариацию. Ограничьте максимальное запаздывание 50% сигнала, чтобы получить хорошую оценку перекрестной ковариации.

maxlags = numel(temp)*0.5;
[xc,lag] = xcov(temp,maxlags);         

[~,df] = findpeaks(xc,'MinPeakDistance',5*2*24);
[~,mf] = findpeaks(xc);

figure
plot(lag/(2*24),xc,'k',...
     lag(df)/(2*24),xc(df),'kv','MarkerFaceColor','r')
grid on
xlim([-15 15])
xlabel('Time (days)')
title('Auto-covariance')

Наблюдайте доминантные и незначительные колебания автоковариации. Доминантные и второстепенные пики кажутся равноудаленными. Чтобы проверить, есть ли они, вычислите и выведите на график разницу между расположениями последующих пиков.

cycle1 = diff(df)/(2*24);
cycle2 = diff(mf)/(2*24);

subplot(2,1,1)
plot(cycle1)
ylabel('Days')
grid on
title('Dominant peak distance')
subplot(2,1,2)
plot(cycle2,'r')
ylabel('Days')
grid on
title('Minor peak distance')

mean(cycle1)

ans = 7

mean(cycle2)

ans = 1

Малые пики указывают 7 циклов в неделю, а доминирующие пики указывают 1 цикл в неделю. Это имеет смысл, учитывая, что данные поступают из контролируемого температурой здания по 7-дневному календарю. Первый 7-дневный цикл показывает, что существует недельное циклическое поведение температуры здания, при котором температура снижается в выходные дни и возвращается к норме в течение недельных дней. Поведение 1-дневного цикла указывает на то, что существует и суточное циклическое поведение - температуры ниже ночью и повышаются днем.