% correlation_example.m

%

% This MATLAB script file demonstrates

% the computation of the correlation

% between two signals.

%

% K. Beers

% MIT ChE. 8/15/2002

% First, compute the correlation between

% two periodic functions.

% Generate time sample information

N = 2^15;

P = 2*pi;

t_val = linspace(0,4*pi,N);

dt = t_val(2) - t_val(1);

t0 = t_val(1);

% Generate f(t)

f = zeros(size(t_val));

for k=1:10

    f = f + (1/sqrt(k))*sin(k*t_val);

end

% Generate g(t)

g = zeros(size(t_val));

phase = 0.5;

g = cos(3.*t_val + phase) + ...

    sin(5.*t_val + phase) + cos(12.*t_val);

% Compute correlation using FFT techniques

f_FT = fft(f);

f_PS = real(conj(f_FT).*f_FT);

g_FT = fft(g);

g_PS = real(conj(g_FT).*g_FT);

gCORRf_FT = g_FT.*conj(f_FT);

gCORRf_PS = conj(gCORRf_FT).*gCORRf_FT;

gCORRf = real(ifft(gCORRf_FT));

% Calculate sampled frequencies

omega = zeros(size(t_val));

omega_crit = pi/dt;

P = (t_val(N) - t_val(1))/2;  % half-period of sample

d_omega = pi/P;

omega = linspace(0,2*omega_crit-d_omega,N);

% Show signals and their correlation

figure;

omega_plot_max = 15;

% f

subplot(3,2,1); plot(t_val,f);

xlabel('t'); ylabel('f(t)');

title('Computation of Corr(g,f)');

axis tight;

% g

subplot(3,2,2); plot(t_val,g);

xlabel('t'); ylabel('g(t)');

axis tight;

% f_PS

subplot(3,2,3);

plot(omega,f_PS);

axis([0 omega_plot_max min(f_PS) max(f_PS)]);

xlabel('\omega'); ylabel('|F(\omega)|^2');

% g_PS

subplot(3,2,4);

plot(omega,g_PS);

axis([0 omega_plot_max min(g_PS) max(g_PS)]);

xlabel('\omega'); ylabel('|G(\omega)|^2');

% gCORRf_PS

subplot(3,2,5);

plot(omega,gCORRf_PS);

axis([0 omega_plot_max min(gCORRf_PS) max(gCORRf_PS)]);

xlabel('\omega'); ylabel('|G F^*|^2');

% gCOORf

subplot(3,2,6);

plot(t_val,gCORRf);

xlabel('t'); ylabel('Corr(g,f)');

axis tight;

% Now, compute convolution of g and f, and compare

% to correlation.

gCONVf_FT = g_FT.*f_FT;

gCONVf_PS = conj(gCONVf_FT).*gCONVf_FT;

gCONVf = real(ifft(gCONVf_PS));

figure;

% gCORRf_PS

subplot(2,2,1);

plot(omega,gCORRf_PS);

axis([0 omega_plot_max min(gCORRf_PS) max(gCORRf_PS)]);

xlabel('\omega'); ylabel('|G F^*|^2');

% gCOORf

subplot(2,2,2);

plot(t_val,gCORRf);

xlabel('t'); ylabel('Corr(g,f)');

axis tight;

% gCONVf_PS

subplot(2,2,3);

plot(omega,gCONVf_PS);

axis([0 omega_plot_max min(gCONVf_PS) max(gCONVf_PS)]);

xlabel('\omega'); ylabel('|G F|^2');

% gCONVf

subplot(2,2,4);

plot(t_val,gCONVf);

xlabel('t'); ylabel('g*f');

axis tight;

% place title manually

gtext('Comparison of correlation and convolution');