MATLAB学习help之——Measuring Signal Similarities from Signal Processing Toolbox

这个例子是用来查看两个信号的相似性，可以针对不同的长度，不同的采样频率的两个信号。

Step1.
加载两个不同采样频率的数据。

% Load data
load relatedsig.mat;

figure
ax(1) = subplot(311);
plot((0:numel(T1)-1)/Fs1,T1,'k');
ylabel('Template 1');
grid on
ax(2) = subplot(312);
plot((0:numel(T2)-1)/Fs2,T2,'r');
ylabel('Template 2');
grid on
ax(3) = subplot(313);
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])

信号如下

图片.png

信号1和信号2为模板信号，信号3为带搜寻的信号，三个信号的频率如下

图片.png

对于不同长度的信号无法执行差分，但是可以通过提取公共部分实现。而互相关则对长度没有要求，但需要保证是要有一致的采样频率。所以可以重采样到一个低频率上， 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

Step2.
寻找信号。
把信号S分别和信号T1，T2取互相关。

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

figure
ax(1) = subplot(211);
plot(lag1/Fs,C1,'k');
ylabel('Amplitude');
grid on
title('Cross-correlation between Template 1 and Signal')
ax(2) = subplot(212);
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 ])

如下

图片.png

从图上看出，目标信号和T1信号弱相关，但可在T2信号中找到，执行如下指令

图片.png

可知信号峰值在61ms

Step3.
测量信号间的延迟，并对齐

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

图片.png

对S1和S2，S1和S3取互相关，则互相关的最大值的点可以表明时间的延迟或滞后

[C21,lag1] = xcorr(s2,s1);
[C31,lag2] = xcorr(s3,s1);

figure
subplot(211);
plot(lag1,C21/max(C21));
ylabel('C21');
grid on
title('Cross-Correlations')
subplot(212);
plot(lag2,C31/max(C31));
ylabel('C31');
grid on
xlabel('Samples')

[~,I1] = max(abs(C21));     % Find the index of the highest peak
[~,I2] = max(abs(C31));     % Find the index of the highest peak
t21 = lag1(I1)              % Time difference between the signals s2,s1
t31 = lag2(I2)              % Time difference between the signals s3,s1

图片.png

t21表明S2滞后S1 350ms， t31表明 S3超前S1 150ms，根据计算结果对其信号如下

s2 = [zeros(abs(t21),1);s2];
s3 = s3(t31:end);

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

图片.png

Step4.
比较信号的频域

Fs = FsSig;         % Sampling Rate

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

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

图片.png

mscohere 函数可以用来比较两个信号的谱特性。

[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(211);
plot(f,Cxy);
title('Coherence Estimate');
grid on;
hgca = gca;
hgca.XTick = f(locs);
hgca.YTick = .75;
axis([0 200 0 1])
subplot(212);
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])

图片.png

可以看到S1，S2信号在35和165Hz附件有强相关内容
在35hz ，相位落后90，在165Hz，落后90度

Step5.
寻找信号的周期性

load officetemp.mat  % Load Temperature Data
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

图片.png

每30分钟记录一次数据，大约16.5周

xcov 在计算互相关的时候，会先去除信号的平均值，返回的是互协方差。

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

图片.png

然后看信号峰值的周期性

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

subplot(211);
plot(cycle1);
ylabel('Days');
grid on
title('Dominant peak distance')
subplot(212);
plot(cycle2,'r');
ylabel('Days');
grid on
title('Minor peak distance')

mean(cycle1)
mean(cycle2)

图片.png

得到的周期时间与数据本身是一致的

总结：

学会几个函数的用法
rat，用于有理数近似
resample，重采样
xcorr，取互相关
periodogram，谱密度分析函数
mscohere，谱相干性分析
cpsd，交叉谱密度分析
angle，取相位
findpeaks，查找峰值点
xcov，取互协方差
diff，差分
不同采样信号频率的相关分析，重采样，数据平移对齐操作
谱密度相关的分析
信号周期性分析。这个地方为什么要用到互协方差暂时没有搞明白。