function [f,pu] = pwrspect(t,u);
% function [f,pu]=pwrspect(t,u);
% input is u(t), output is pu(f)
% pu(f) is: power spectrum of input signal u(t)
% t,u,f, and pu are vectors.  Function will also dispay peak freq.
n0 = length(t);
if rem(n0,2) == 1;
 t(n0) = [];
 u(n0) = [];
end;
av=mean(u);
n = length(t);
nh = n/2;
dt = (t(n) - t(1))/n;
f0 = .5/dt;
f = f0*(0:nh)/nh;
trans = fft(u-av);
pu = 2*trans.*conj(trans)/n^2;
pu(nh+2:n) = [];
pu(2:nh) =  2*pu(2:nh);
%print peak freq
[pumax,m]=max(pu);
disp('peak frequency (Hz) = ')
disp(f(m))
%plot results
plot(f,pu)
title('Power Spectral Density');
ylabel('Normalized Power');
xlabel('Frequency in Hz');