How can I get the residuals of phase between a wave and a modelled wave?

I'm using the following script to build a wave and with the FFT extract the frequency information to build a wave based on that frequency. I managed to make it so that I can "predict" 250 ms of the signal based on the previous 250 ms. However, to validate my model I need to compare the phase and frequency (it's the most important part of my project) of the prediction with the actual data. I have, however no idea of how to accomplish this.

The script below relies on the script "phasereset", which is used to simulate EEG data. The script "phasereset" depends on the eeglab package.

The phasereset script is available at: http://www.cs.bris.ac.uk/~rafal/phasereset/

% create time vector t = 1/250:1/250:1000*(1/250);

%simulate EEG signal

mynoise=noise(1000,1,250);

alpha = phasereset (1000, 1, 250, 8, 13, 0);

beta = phasereset (1000, 1, 250, 13, 30, 0);

gamma = phasereset (1000, 1, 250, 30, 80, 0);

delta = phasereset (1000, 1, 250, .1, 4, 0);

theta = phasereset (1000, 1, 250, 4, 8, 0);

myEEG = alpha+beta+gamma+delta+theta+mynoise;

% define base epoch, test epoch and epoch to mark alpha phases

baseEpoch = myEEG(1:250); % base epoch: we base the alpha model on this

testEpoch = myEEG(251:500); % we validate the alpha model on this dataset

markEpoch = myEEG(501:750); % ...and project our markers on ths dataset

% plot the epochs

h = figure; subplot(2,3,1); hold on;

plot(t(1:250),baseEpoch); title('Baseline epoch');

subplot(2,3,2);

plot(t(251:500),testEpoch); title('Test epoch');

subplot(2,3,3);

plot(t(501:750),markEpoch); title('Epoch to mark');

% now we need to construct our alpha model

%

% a model sine wave equation is:

%

% y(t) = A * sin(2 * pi * f * t + phi)

%

% where A = amplitude, f = frequency, and phi = phase

%

% the parameters can be extracted from the base epoch

% do FFT on the base epoch

Fs = 250; % sample rate in Hz

L = length(baseEpoch); % length of the base epoch in samples

NFFT = 2^nextpow2(L);

Y = fft(testEpoch,NFFT)/L;

f = Fs/2*linspace(0,1,NFFT/2+1); % f is a vector with all frequencies on which we now have information

i = find(f>8 & f<13); % get the indices to the alpha band

% Y is a vector of complex numbers: every value Y(i) is a polar coordinate

% giving you the amplitude of the frequency f(i) and the phase of frequency

% f(i) at the onset of the signal. The amplitude is given by the complex

% modulus or magnitude of Y(i) and phase is given by the angle of Y(i)

% first, get the maximum amplitude in the alpha band and the associated % frequency

[dum, peakAlpha_i] = max(2*abs(Y(i)));

% peakAlpha_i now points to the largest amplitude in the alpha band.

% However, f is defined over the entire frequency spectrum, so we need to

% add an offset to get the actual index number (which is (i(1) - 1))

mod_frequency = f(peakAlpha_i + (i(1)-1)); % frequency of our alpha model

% we also need an estimate of amplitude

mod_amplitude = 2*abs(Y(peakAlpha_i + (i(1)-1))); % amplitude component for alpha model

% and finally phase

mod_phase = angle(Y(peakAlpha_i +(i(1)-1))) + (pi); % phase of our alpha model (in rad)

% plot the alpha band

subplot(2,3,4); hold on;

plot(f,2*abs(Y(1:length(f))));

xlim([8 13]);

line([mod_frequency mod_frequency],[0 1]);

ylim([0 1]);

title('FFT, zoomed in on alpha band');

% now we can build our alpha model:

test_time = t(251:500);

test_alpha = mod_amplitude * sin(2*pi*mod_frequency*test_time + mod_phase);

subplot(2,3,5); hold on;

plot(test_time, alpha(251:500));

title('Test alpha and forward model');

plot(test_time, test_alpha, 'k');