ISU EE 573 Matlab Programs S98

ISU STAT 333 SPRING 1998 Matlab Programs

MV/AR program

% Program Name:  SC333.m
% =====================================================
% REQUIRED RESIDENT INPUT: 
% vector 'ts' containing data
% 'fsamp' = sample frequency associated with ts file
% 'ordsve'=[n1,n2,n3,n4] MV(n) orders

% EXAMPLE:
% ordsve=[5 10 20 40];
%
%   ===================================================
% OTHER REQUIRED INPUT PROVIDED BY THIS PROGRAM:
maxorder=max(ordsve); 
% max order to be computed
nfft=2048; %number of FFT points for desired freq. resolution
%  ***************************************************************
[r,c]=size(ts);   % make sure ts is a column vector
if r == 1
  ts=ts';
end
ts=ts-mean(ts);   % remove mean of ts
%************************************************************
% COMPUTE AUTOCORRELATION SEQUENCE
% ================================   
N=length(ts);   
r(1)=(ts' * ts)/N;   
xvec=[ts ; zeros(N,1)];   
for t=1:maxorder      
y=[zeros(t,1) ; 
xvec(1:(2*N-t))];      
r(t+1)=(xvec' * y) / (N-t);   
end   
%***********************************************************  
%  COMPUTE LEVINSON POLYNOMIALS
%  AS A FUNCTION OF exp(iw)
%   ===================================   
E2=[];   
E2(1)=r(1);   
Alpha=1.0;   
Aall=[Alpha; zeros(maxorder,1)];
for n=1:maxorder      
  rflip=fliplr(r(1:n+1));      
  Alpha=[Alpha; 0.0];      
  del=rflip * Alpha;      
  Alpha=Alpha - (del/E2(n)) * flipud(Alpha);      
  E2(n+1) = E2(n) - (del^2)/E2(n);      
  Aall=[Aall,[Alpha; zeros(maxorder-n,1)]];   
end
Aall=flipud(Aall);
%*******************************************************
%  COMPUTE MOD-SQD. FFTS OF POLYNOMIALS
%  ===================================
sumphi2=zeros(nfft,1);
for order=0:maxorder   
 order   
 phi=(fft(Aall(:,order+1),nfft)) ./ sqrt(E2(order+1));   
 phimod2=(abs(phi)) .^ 2;   
 sumphi2=sumphi2 + phimod2;   
%****************************
% SEARCH FOR SUMS TO BE SAVED
  if order == ordsve(1)      
    mv1= 1 ./ sumphi2;      
    ar1= 1 ./ phimod2;   
  elseif order == ordsve(2)      
    mv2= 1 ./ sumphi2;      
    ar2= 1 ./ phimod2;   
  elseif order == ordsve(3)      
    mv3= 1 ./ sumphi2;      
    ar3= 1 ./ phimod2;   
  elseif order == ordsve(4)      
    mv4= 1 ./ sumphi2;      
    ar4= 1 ./ phimod2;   
  end
end
%****************************
NAr=[ar1 ,ar2, ar3, ar4];
NAr2=NAr(1:nfft/2,:);
NMv=[mv1 ,mv2, mv3, mv4];
NMv2=NMv(1:nfft/2,:);
faxis=0:fsamp/nfft:(fsamp/2)-(fsamp/nfft);
figure(1)
plot(faxis,10*log10(NAr2))
title('AR SPECTRA')
xlabel('FREQUENCY');
ylabel(' dB')
pause
figure(2)
plot(faxis,10*log10(NMv2))
title('MV SPECTRA')
xlabel('FREQUENCY');
ylabel(' dB')