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

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

% Load data
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 функция применяет сглаживающийся (lowpass) КИХ-фильтр к сигналу во время процесса передискретизации.

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

Наблюдайте доминирующие и незначительные колебания автоковариации. Доминирующий и незначительный 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-дневное поведение цикла указывает, что существует также ежедневное циклическое поведение - температуры ниже в течение ночи и увеличения в течение дня.