MATLAB Answers

Tiera-Brandy

changing time scale in cross correlation

Asked by Tiera-Brandy
on 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 Comment

Daniel
on 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.

Products

No products are associated with this question.

0 Answers

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply today