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

图片.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

图片.png

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

总结:

  1. 学会几个函数的用法
    rat,用于有理数近似
    resample,重采样
    xcorr,取互相关
    periodogram,谱密度分析函数
    mscohere,谱相干性分析
    cpsd,交叉谱密度分析
    angle, 取相位
    findpeaks, 查找峰值点
    xcov,取互协方差
    diff, 差分

  2. 不同采样信号频率的相关分析,重采样,数据平移对齐操作

  3. 谱密度相关的分析

  4. 信号周期性分析。这个地方为什么要用到互协方差暂时没有搞明白。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,688评论 1 330
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,559评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,749评论 0 226
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,581评论 0 191
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,741评论 3 271
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,684评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,122评论 2 292
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,847评论 0 182
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,441评论 0 228
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,939评论 2 232
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,333评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,783评论 2 236
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,275评论 3 220
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,830评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,444评论 0 180
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,553评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,618评论 2 249