% freqanal_paramtreadmill.m
% William C Rose
% 2/5/2006,3/21/2006,5/14/2007,6/5/2007,20080131,20080402,20080415
% Requires Signal Processing Toolbox and external function mdl(m) by WCR.
% Read text file with x(k)=column of numbers data sampled at uniform rate.
% Does parametric spectrum analysis (AR model) of x(k), using forward-backward 
% algorithm (which I've shown to be same as modified covariance method).  
% Tries many orders, finds orders with best AIC, best MDL.  Used to do FPE.
% Best FPE order = best AIC order in all data I've tried including motor
% units, heart rate, RR intervals.  Now FPE is gone.
% Evaluates whiteness of residuals. Plot residual autocorrelation.
% Edited 20080415: adjust frequencies for use wil treadmill force data from
% grf_z_treadmill_1khz_raw.txt.

clear;
infilespec = input('Enter name of text file containing signal: ','s');
[pathstr,name,ext,versn] = fileparts(infilespec); %use later in plot title
y1=load(infilespec);
N=length(y1);
fprintf('Number of points = %d.\n',N);
Fs = input('Enter signal sampling rate, in Hertz: ');
% Remove linear trend from the full length data record.
t = (0:1/Fs:(N-1)/Fs)';     % make a vector of times
p = polyfit(t,y1,1);        % find coeffs of straight line thru the data
y = y1-polyval(p,t);        % remove linear trend from the data
Data = iddata(y,[],1/Fs);   % set up data for modeling

% Loop to find AIC, MDL at many different orders, coarsely first
p=[6:2:28 32:4:200];
aic_a=p; mdl_a=p;
for i = 1:length(p)
    m_fb=ar(Data,p(i),'fb');
    aic_a(i)=aic(m_fb);
    mdl_a(i)=mdl(m_fb);
end;
[aicbest,i_aic]=min(aic_a);       % model with best AIC so far
[mdlbest,i_mdl]=min(mdl_a);       % model with best MDL so far
% loop again now fine search
p=[p(i_aic)-4:p(i_aic)+4 p(i_mdl)-4:p(i_mdl)+4];
aic_a=p; mdl_a=p;
for i = 1:length(p)
    m_fb=ar(Data,p(i),'fb');
    aic_a(i)=aic(m_fb);
    mdl_a(i)=mdl(m_fb);
end;
[aicbest,i_aic]=min(aic_a);       % model with best AIC
[mdlbest,i_mdl]=min(mdl_a);       % model with best MDL
fprintf('Best AIC: order %d, AIC=%8.4g\n',p(i_aic),aicbest);
fprintf('Best MDL: order %d, MDL=%8.4g\n',p(i_mdl),mdlbest);

% Plot spectra using best orders
f=[Fs/N:Fs/N:Fs/2];
[Pxxaic f]=pmcov(y,p(i_aic),f,Fs);      % order with best AIC
figure;
subplot(3,1,1); plot(f,Pxxaic,'-b.');
titlestring=sprintf('freqanal_param: %s%s (Fs=%.0f Hz, N=%d)  Best AIC order=%d',name,ext,Fs,N,p(i_aic));
title(titlestring,'Interpreter','none'); 
ylabel('POWER SPECTRAL DENSITY');
grid on; % xlim([0 5]);

[Pxxmdl f]=pmcov(y,p(i_mdl),f,Fs);      % order with best MDL
subplot(3,1,2);  plot(f,Pxxmdl,'-b.');
titlestring=sprintf('Best MDL order=%d',p(i_mdl));
title(titlestring,'Interpreter','none'); 
ylabel('POWER SPECTRAL DENSITY');
grid on; % xlim([0 5]); 

[Pxx120 f]=pmcov(y,120,f,Fs);           % order=120
subplot(3,1,3);  
plot(f,Pxx120,'-b.');
titlestring=sprintf('Order=120');
title(titlestring,'Interpreter','none'); 
ylabel('POWER SPECTRAL DENSITY'); xlabel('FREQUENCY (Hz)');
grid on; % xlim([0 5]);

% Plot overlaid spectra using different orders
figure;
for i=1:9                               % was i=1:5 for EB2008, 20080402
    pp=(i+1)*20;                        % orders 40-200, was pp=40*i for EB2008
    [Pxx(:,i) f]=pmcov(y,pp,f,Fs);
    legstr{i}=sprintf('Order=%d',pp);     % cell array of strings for legend
end;
plot(f,Pxx(:,1),f,Pxx(:,2),f,Pxx(:,3),f,Pxx(:,4),f,Pxx(:,5),f,Pxx(:,6),f,Pxx(:,7),f,Pxx(:,8),f,Pxx(:,9));
%plot(f,Pxx(:,1),f,Pxx(:,2),f,Pxx(:,3),f,Pxx(:,4),f,Pxx(:,5));
titlestring=sprintf('freqanal_param: %s%s (Fs=%.0f Hz, N=%d)',name,ext,Fs,N);
title(titlestring,'Interpreter','none'); 
ylabel('POWER SPECTRAL DENSITY'); xlabel('FREQUENCY (Hz)');
legend(legstr);
grid on; xlim([0 5]);
mrkr = 'osx.+*dv^<>ph';             % array of marker styles
fc = get(gcf,'children');           % handles to legend & axis
ac = flipud(get(fc(2),'children')); % handles to lines in plot, reordered
for i = 0:length(ac)-1
    set(ac(i+1),'marker',mrkr(mod(i,length(mrkr))+1))
end
%setprinttemplate(gcf,readDefProps('portrait bothalf color'));
%print       % causes computer crash when print to Dell 926 AIO printer

% Prediction error whiteness test (PEWT)
% If the model residuals are white (i.e. uncorrelated, i.e. autocorr=0 
% at all lags) then about 1/100 should be > 99% C.I. for autocorr=0 and
% 5/100 should be > 95% CI for autocorr=0. 
% Plot autocorr for lags 1-100 and 99% CI for autocorr=0.
pp=120;
m=ar(Data,pp,'fb');
N=m.es.DataLength;          %N=number of data points fitted
[e1 r1]=resid(m,Data,101);  %first 100 lags of autocorrelation func
thresh95=1.96/sqrt(N);
thresh99=2.58/sqrt(N);
Nover95=length(find(abs(r1(2:101))/r1(1)>thresh95));
Nover99=length(find(abs(r1(2:101))/r1(1)>thresh99));
figure;
resid(m,Data,101);          %first 100 lags of autocorrelation func
tstr1=sprintf('Residual Autocorrelation, order %d: %s',pp,infilespec);
title(tstr1,'Interpreter','none');
ylim([-3*thresh99 3*thresh99]);
legend(['99% CI'],['autocorr: ',num2str(Nover99),'/100 outside 99% CI']);
fprintf('Model order %d: %d/100 outside 95%% CI. ',pp,Nover95);
if Nover95>10
    fprintf('Probably not white.\n') 
else
    fprintf('Probably white.\n') 
end;
fprintf('Model order=%d: %d/100 outside 99%% CI. ',pp,Nover99);
if Nover99>4 
    fprintf('Probably not white.\n')
else
    fprintf('Probably white.\n')
end;