function spsd(A,Bw,C,Sw,Sv,omega1,omega2)
% SPSD  System Power Spectral Density Plot.
% spsd(A,Bw,C,Sw,Sv,omega1,omega2) plots psd of all outputs on a single
% graph.  The frequency axis is from omega1 to omega2 rads per sec
% Plant model is 
%  .
%  x = Ax + Bw w
%  y = Cx + v
%
% R.J. Vaccaro 4/09
%
d1=omega1;d2=omega2;

npoints=round((d2-d1)*200);
w=linspace(d1,d2,npoints);
z=j*w;

% Evaluate L(z) over the frequencies in z vector.

nz = length(z);

[t,a] = balance(A);
b = t \ Bw;
c = C * t;

[p,a] = hess(a);
b = p' * b;
c = c * p;
[nop,xx]=size(C);
[n,m]=size(b);
L=zeros(nz,nop);
for kkk=1:m
    g = ltifr(a,b(:,kkk),z);
    for kk=1:nop
        for k=1:nz
            L(k,kk)=L(k,kk)+real(c(kk,:)*g(:,k)*Sw(kkk,kkk)*.....
                (g(:,k)'*c(kk,:)'))+Sv(kk,kk);
        end
    end
end
hold off
plot(w,L);
grid
title('Power Spectral Density');
xlabel('Frequency (rad/sec)');
