%
%  The following code segment is used to estimate the psd three ways.
%  Note that we do not have an estimate of the variance for the ARMA
%  method, though we could get it with more work....
%
%  Now use the data to estimate the parameters of the ARMA model
%
%  put the data in the correct structure
%
  y = iddata(x',[], Ts);
%
%  guess the number of AR and MA parts
%
  na = 4; % number of coefficients for the AR part (how many a's)
  nc = 1; % number of coefficients for the MA part, minus 1 (how many b's --ignoring b0)
%
%  determine the model parameters
%
  sys = armax(y,[na nc])
  [A,B,C,D] = polydata(sys);
  getcov(sys)
%
%  make this look the way we want
%
  B = C;
%
%  Now get ready to plot
%
  figure;
%
%  first the ARMA based results
%
  [H_ARMA,F_ARMA] = freqz(B,A,1000,1/Ts);
%
%  Now use the FFT Welch method
%
  [Pxx,F_Welch] = pwelch(x,[],[],[],1/Ts);
%
%  Now computer the FFT of the autocorrelation
%  
  Sxx = (1/r)*abs(fft(Rxx));
  half_Sxx = Sxx(1:N/2);
  F_Sxx = linspace(0,0.5/Ts,length(half_Sxx));
%
%  YOu might want to scale
%
  scale_P = mean(H_ARMA)/mean(Pxx);
  scale_S = mean(H_ARMA)/mean(half_Sxx);
%   scale_W = 1.0;
%   scale_S = 1.0;
%  Now plot the results
%
  hold on
  plot(F_Welch,20*log10(Pxx*scale_P),'r-','LineWidth',1.5);
  plot(F_ARMA,20*log10(H_ARMA),'b--','LineWidth',1.0);
  plot(F_Sxx,20*log10(half_Sxx*scale_S),'m-','LineWidth',1.0);
  hold off
  legend('Welch (FFT)','ARMA','Sxx');
  grid; xlabel('Frequency (Hz)'); ylabel('PSD (dB)');
%

