www.gusucode.com > signal 案例源码程序 matlab代码 > signal/SignalSimilaritiesExample.m
%% Measuring Signal Similarities % This example shows how to measure signal similarities. It will help you % answer questions such as: How do I compare signals with different lengths % or different sampling rates? How do I find if there is a signal or just % noise in a measurement? Are two signals related? How to measure a delay % between two signals (and how do I align them)? How do I compare the % frequency content of two signals? Similarities can also be found in % different sections of a signal to determine if a signal is periodic. % Copyright 2012-2014 The MathWorks, Inc. %% Comparing Signals with Different Sampling Rates % Consider a database of audio signals and a pattern matching application % where you need to identify a song as it is playing. Data is commonly % stored at a low sampling rate to occupy less memory. % 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]) %% % The first and the second subplot show the template signals from the % database. The third subplot shows the signal which we want to search for % in our database. Just by looking at the time series, the signal does not % seem to match to any of the two templates. A closer inspection reveals % that the signals actually have different lengths and sampling rates. [Fs1 Fs2 Fs] %% % Different lengths prevent you from calculating the difference between two % signals but this can easily be remedied by extracting the common part of % signals. Furthermore, it is not always necessary to equalize lengths. % Cross-correlation can be performed between signals with different % lengths, but it is essential to ensure that they have identical sampling % rates. The safest way to do this is to resample the signal with a lower % sampling rate. The |resample| function applies an anti-aliasing(low-pass) % FIR filter to the signal during the resampling process. [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 %% Finding a Signal in a Measurement % We can now cross-correlate signal S to templates T1 and T2 with the % |xcorr| function to determine if there is a match. [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 ]) %% % The first subplot indicates that the signal and template 1 are less % correlated while the high peak in the second subplot indicates that % signal is present in the second template. [~,I] = max(abs(C2)); SampleDiff = lag2(I) timeDiff = SampleDiff/Fs %% % The peak of the cross correlation implies that the signal is present in % template T2 starting after 61 ms. In other words, signal T2 leads signal % S by 499 samples as indicated by SampleDiff. This information can be used % to align the signals. %% Measuring Delay Between Signals and Aligning Them % Consider a situation where you are collecting data from different % sensors, recording vibrations caused by cars on both sides of a bridge. % When you analyze the signals, you may need to align them. Assume you have % 3 sensors working at same sampling rates and they are measuring signals % caused by the same event. 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') %% % We can also use the |finddelay| function to find the delay between two % signals. t21 = finddelay(s1,s2) t31 = finddelay(s1,s3) %% % t21 indicates that s2 lags s1 by 350 samples, and t31 indicates that s3 % leads s1 by 150 samples. This information can now used to align the 3 % signals by time shifting the signals. We can also use the |alignsignals| % function directly to align the signals which aligns two signals by % delaying the earliest signal. 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') %% Comparing the Frequency Content of Signals % A power spectrum displays the power present in each frequency. Spectral % coherence identifies frequency-domain correlation between signals. % Coherence values tending towards 0 indicate that the corresponding % frequency components are uncorrelated while values tending towards 1 % indicate that the corresponding frequency components are correlated. % Consider two signals and their respective power spectra. 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| function calculates the spectral coherence between the two % signals. It confirms that sig1 and sig2 have two correlated components % around 35 Hz and 165 Hz. In frequencies where spectral coherence is high, % the relative phase between the correlated components can be estimated % with the cross-spectrum phase. [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]) %% % The phase lag between the 35 Hz components is close to -90 degrees, and % the phase lag between the 165 Hz components is close to -60 degrees. %% Finding Periodicities in a Signal % Consider a set of temperature measurements in an office building during % the winter season. Measurements were taken every 30 minutes for about % 16.5 weeks. % 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 %% % With the temperatures in the low 70s, you need to remove the mean to % analyze small fluctuations in the signal. The |xcov| function removes the % mean of the signal before computing the cross-correlation. It returns the % cross-covariance. Limit the maximum lag to 50% of the signal to get a % good estimate of the cross-covariance. 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') %% % Observe dominant and minor fluctuations in the auto-covariance. Dominant % and minor peaks appear equidistant. To verify if they are, compute and % plot the difference between the locations of subsequent 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) mean(cycle2) %% % The minor peaks indicate 7 cycles/week and the dominant peaks indicate 1 % cycle per week. This makes sense given that the data comes from a % temperature-controlled building on a 7 day calendar. The first 7-day % cycle indicates that there is a weekly cyclic behavior of the building % temperature where temperatures lower during the weekends and go back to % normal during the week days. The 1-day cycle behavior indicates that % there is also a daily cyclic behavior - temperatures lower during the % night and increase during the day.