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

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

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

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

[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)')

The 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-х нужно удалить среднее для анализа небольших колебаний сигнала. The 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')

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

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

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