% Example of 1D parametric spectral analysis using autoregressive process.
%
% Uses the first row of the first frame of video sequence as example
% signal, estimates a pth order autoregressive process for the data (using
% covariance method), and converts this to a magnitude spectrum.  The model
% can also be used to synthesise new data samples with the same correlation
% structure.  
%
% Notes:
% 1. The mean is removed from the data before the estimation is done.  Is
%    there any good motivation for this?
% 2. The methods demonstrated in this code extend quite simply to 2D, in
%    which case the function im2col is still appropriate.  
% 3. When extended to 3D, namely 2D + time, the model is sometimes called a
%    spatio-temporal autoregressive (STAR) model.  The function im2colp
%    extends the standard function im2col to the higher-dimensional case.

% Get data
%mov = aviread('sea_waves.avi');
%im = mov(1).cdata;  % first frame
%mov = VideoReader('sea_waves.avi');
%im = readFrame(mov)
im = imread('sea_waves2.png');
dat1 = double(im(1,:));  % first row
dat1 = dat1 - mean(dat1);  % ???
fh = figure;  plot(dat1);  title('Actual data (for training)');

% Get samples for linear prediction.  To understand the function im2dat
% used below, look carefully at the inputs and outputs of the command
% im2col(1:10,[1 5])
p = 4;  % model order (effect?  choice?)
X0 = flipud(im2col(dat1,[1 p+1]));  % data matrix
X = X0(2:end,:)';  y = X0(1,:)';  % observation vectors and target values

% Solve for model
a = X\y;  % equivalent to a = inv(X'*X)*X'*y
w = y-X*a;  % error sequence (should be white if model is good)
sigsq = var(w);  % sample variance of white noise driving sequence
fh = figure;  plot(w);  title('Estimated white noise sequence from training data');

% Spectrum estimate (is the DC gain correct?)
wv = linspace(0,2*pi,200);  % frequencies to calculate
denv = ones(size(wv));
for k=1:p
  denv = denv - a(k)*exp(-1j*wv*k);
end
Hv = sqrt(sigsq)./denv;  % transfer function values (complex)
Hv2 = freqz(sqrt(sigsq),[1 -a'],wv);  % should be same as Hv above
fh = figure;  
subplot(2,1,1);  plot(wv,abs(Hv));  axis tight;  xlabel('\omega');  ylabel('|H(e^{j \omega})|');  
subplot(2,1,2);  plot(wv,angle(Hv));  axis tight;  xlabel('\omega');  ylabel('\angle H(e^{j \omega})');

% Synthetic sample from estimated model
wnsig = sqrt(sigsq)*randn(size(dat1));  % new white noise driving sequence
aa = [1 -a'];  bb = 1;  % compare specified model to expression in "doc filter"
ssig = filter(bb,aa,wnsig);  % note that inverse is wnsig = filter(aa,bb,ssig)
fh = figure;  plot(ssig);  title('Synthetic data from model');