changing time scale in cross correlation

8 Ansichten (letzte 30 Tage)
Tiera-Brandy
Tiera-Brandy am 9 Jan. 2012
Hi all, I'm running a 1D cross correlation on two timeseries that are given in years, however I want to show a lag in months not years. I'm not sure where to change the script to do this. The following is my script to create the correlation and lag:
t1m=time(10:341); %Excludes ends of time length to
s1m=ssha_area(10:341); %better fit data line
t2m=ctime(10:160);
s2m=chlo_area(10:160);
N1=length(t1m); % sample size
N2=length(t2m); %
N21=ceil(N1/2); % half the number of obs(rounded)
N22=ceil(N2/2); %
ts1 = t1m(1); %
ts2 = t2m(1); %
tlast1 = t1m(end); % last data point
tlast2 = t2m(end); %
ti1 = (tlast1-ts1)/(N1-1); % mean sampling interval
ti2 = (tlast2-ts2)/(N2-1); %
te1 = tlast1+ti1; % end plus one increment
te2 = tlast2+ti2; %
%
% frequency vector
%
fmin1=0; % the min freq is always 0
fmin2=0; %
fi1=1/(te1-ts1); % the lowest non-zero freq in data
% is always 1 divided by the length
% of the data set
fi2=1/(te2-ts2); %
fs1=1/ti1; % the sampling frequency
fs2=1/ti2; %
fmax1=fs1-fi1; % before applying the fftshift, the
% maximum frequency is always the
% sampling frequency minus the
fmax2=fs2-fi2; % frequency increment
fpre1=(fmin1:fi1:fmax1)'; % frequency grid before fftshift
fpre2=(fmin2:fi2:fmax2)'; %
f1=[fpre1(N21+1:N1)-fs1;fpre1(1:N21)];
f2=[fpre2(N22+1:N2)-fs2;fpre2(1:N22)];
% Now apply the Fourier transform
%
S1=fft(s1m,N1);
S2=fft(s2m,N2);
%
% Then apply fftshift to S1 and S2:
%
S1shift=fftshift(S1);
S2shift=fftshift(S2);
%
% and compute the magnitude
%
S1mag=fftshift(abs(S1));
S2mag=fftshift(abs(S2));
% Now plot the original data.
%
figure(1), clf
subplot(2,1,1), grid on, box on, hold on
plot(t1m,s1m,'b--.')
legend('signal 1 unfiltered')
xlabel('time (years)'), ylabel('signal (m)')
subplot(2,1,2), grid on, box on, hold on
plot(f1,S1mag,'b--.')
legend('signal 1 unfiltered signal')
xlabel('frequency (years^{-1})'), ylabel('magnitude (m)')
figure(2), clf
subplot(2,1,1), grid on, box on, hold on
plot(t2m,s2m,'b--.')
legend('signal 2 unfiltered')
xlabel('time (years)'), ylabel('signal (m)')
subplot(2,1,2), grid on, box on, hold on
plot(f2,S2mag,'b--.')
legend('signal 2 unfiltered')
xlabel('frequency (hrs^{-1})'), ylabel('amplitude (m)')
end
This is where the cross correlation begins really
% Cross correlate the two time-series, and find the dominant frequency at which they are correlated, and the phase difference at that frequency.
s1m = s1m-mean(s1m);
s2m = s2m-mean(s2m);
s1m=s1m(1:N2);
[x12pre,lag12pre]=xcorr(s1m,s2m)
Ok I think this is where I would change the script but that would seem to be dependent on my sampling interval(ti1) as well so I'm not sure if there are multiple parts I have to change.
% Convert lag12pre into values of unit years.
%
lag12=((lag12pre*ti1));
%
% Normalise the cross-correlation to fall between -1 and 1
%
[a11,lag11]=xcorr(s1m);
[a22,lag22]=xcorr(s2m);
lagind11=find(lag11==0);
lagind22=find(lag22==0);
x12=x12pre/sqrt(a11(lagind11).*a22(lagind22));
%
% Calculate the cross-spectrum by Fourier transforming the
% cross-correlation.
%
% (ifftshift).
%
x12shift=ifftshift(x12);
lag12shift=lag12-min(lag12);
%
% Create a new frequency vector, corresponding to the
% time-lag axis of the cross-correlation of the signals.
%
N12=length(lag12); % number of observations
NT12=lag12(end)-lag12(1); % length of record (years)
dT12=NT12/(N12-1); % mean sample space (years).
fs12=1/dT12; % sampling frequency (max frequency)
df12=(1/NT12)*(N12-1)/N12; % minimum non-zero frequency.
%
f12pre=0:df12:fs12-df12; % the frequency vector
f12=fftshift(f12pre); % and the vector that corresponds
% to the fftshifted signal
f12(f12>(fs12-df12/2)/2)=f12(f12>(fs12-df12/2)/2)-fs12;
%
% Calculate the Fourier transform, and fftshift it:
%
X12=fft(x12shift);
X12shift=fftshift(X12);
%
% Normalise the magnitude by the length of the record.
%
X12mag=abs(X12shift);
%
% Calculate the angle
X12phs=angle(X12shift);
%
if 1
figure(3), clf
subplot(3,1,1), grid on, box on, hold on
title('Lag Between SSHA and Chlorophyll')
plot(lag12,x12,'b--.')
xlabel('lag (years)'), ylabel('correlation (m^2)')
subplot(3,1,2), grid on, box on, hold on
plot(f12,X12mag,'b--.')
xlabel('frequency (years^{-1})'), ylabel('magnitude of xcorr (m^2)')
subplot(3,1,3), grid on, box on, hold on
plot(f12,X12phs,'b--.')
xlabel('frequency (years^{-1})'), ylabel('lag (rads)')
end
  1 Kommentar
Daniel Shub
Daniel Shub am 10 Jan. 2012
You need to work through your problem enough that you can give a much smaller MWE. Don't give us anything that we don't need ... For example, the question doesn't see to deal with plotting, so get rid of all of that.

Melden Sie sich an, um zu kommentieren.

Antworten (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by