% Script to find power spectral density of u and w for HW1, p3b.

  load('rf18L1.txt');          % Load u and w
  u = rf18L1(:,2); 
  w = rf18L1(:,4); 

  rate = 25;                   % Data rate (Hz)
  nfft = 4096;                 % Samples per data window
  ns = length(w);              % Number of samples
  window = hanning(nfft);      % Use Hanning data taper function
  noverlap = nfft/2;           % Overlap between periodogram windows
  dflag = 'none';              % No detrending
  [Pww,f] = psd(w,nfft,ns,window,noverlap,dflag);
  [Puu,f] = psd(u,nfft,ns,window,noverlap,dflag);
  f = f*rate/ns;               % Convert f to units of s^-1
  Pww = Pww/rate;              % Convert Pww and Puu so 
  Puu = Puu/rate;              % integral of Pww df = var(w), similarly for u.