sprfft.m

% sprfft - Program to compute the power spectrum of a  
% coupled mass-spring system.
clear; help sprfft;  % Clear memory and print header
x = input('Enter initial displacement [x1 x2 x3] - ');  
v = [0 0 0];       % Masses are initially at rest
state = [x v];     % Positions and velocities; used by rk4
tau = input('Enter timestep, tau - ');  % (sec)
k_over_m = 1;      % Ratio of spring const. over mass (1/sec^2)
%%%%%%%%%%%%% Main Loop %%%%%%%%%%%%%%
time = 0;
nstep = 256;       % Number of steps in the main loop
nprint = nstep/16; % Number of steps between printing progress
for istep=1:nstep  
  state = rk4(state,time,tau,'sprrk',k_over_m);  % Runge-Kutta
  time = time + tau;    
  xplot(istep,1:3) = state(1:3);   % Record positions
  tplot(istep) = time;
  if( rem(istep,nprint) < 1 )
    fprintf('Finished %g out of %g steps\n',istep,nstep);
  end
end
f(1:nstep) = (0:(nstep-1))/(tau*nstep);      % Frequency
x1 = xplot(:,1);              % Displacement of mass 1
window = 0.5*(1-cos(2*pi*(1:nstep)/nstep));  % Hanning window
x1w = x1 .* window';
x1fft = fft(x1);              % Fourier transform of displacement
spect = x1fft .* conj(x1fft); % Power spectrum of displacement
x1wfft = fft(x1w);            % Fourier transf. (windowed data)
spectw = x1wfft .* conj(x1wfft); % Power spectrum (windowed data)
% Graph the displacements of the masses
ipr = 1:nprint:nstep;
plot(tplot,xplot(:,1),'-',tplot(ipr),xplot(ipr,1),'o',...
     tplot,xplot(:,2),'-.',tplot(ipr),xplot(ipr,2),'+',...
     tplot,xplot(:,3),'--',tplot(ipr),xplot(ipr,3),'*')
title('Displacement of mass 1- o  2- +  3- *');
xlabel('Time'); ylabel('Displacement');
pause;  % Pause between plots; strike any key to continue
semilogy(f(1:(nstep/2)),spect(1:(nstep/2)),'-',...
         f(1:(nstep/2)),spectw(1:(nstep/2)),'--');
title('Power spectrum (dashed is windowed data)');
xlabel('Frequency'); ylabel('Power');