n = 3;

 m = 0.5;            % .5 oz
 sigma = 5.5;        % 5.5 lbs
 
% units in cm -- reported by Scott
 ell0 = 12;          % 12 cm
 ell1 = 19;          % 19 cm
 ell2 = 19;          % 19 cm
 ell3 = 22;          % 22 cm

% units in cm  -- measured the next morning...
% ell0 = 14;
% ell1 = 20.5;
% ell2 = 21.25;
% ell3 = 22;



 a1 = sigma/2*(1/ell0+1/ell1);
 a2 = sigma/2*(1/ell1+1/ell2);
 a3 = sigma/2*(1/ell2+1/ell3);
 b1 = sigma/(2*ell1);
 b2 = sigma/(2*ell2);
 

 C = .5*m*eye(3);
 A = [ a1 -b1   0 ; 
      -b1  a2 -b2 ;
        0 -b2  a3 ]; 

 scale =  (16)*(1/0.45359237) * (4.4482216162605) * (100);
%         (oz/lb)*(lb/kg)            (N/lbf)        (cm/m)
% scale from Hardesty's calculator = 15690

 B = scale*inv(C)*A;

 [V,D] = eig(B);
 rho = sqrt(diag(D));
 freq = rho/(2*pi);
 fprintf('Predicted natural frequencies:\n')
 disp(freq)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% edited from Sean Hardesty's fftwav.m

 filename = 'bead3a.wav';

% plots time and frequency data for wav file filename
% Sean Hardesty (1 Feb 2006)

 [y,fs]=wavread(filename);
 N=2^floor(log(length(y))/log(2));
 if (N < 0.75*length(y))
      N=length(y);
 end

%use the next line if you have an old version of Matlab
  w=fft(y(1:N,:).*kron(sin((0:(N-1))*pi/N)',ones(1,size(y,2)))); % me 
% w=fft(y(1:N,:).*kron(window(@hamming,N),ones(1,size(y,2))));

 range=(1:(N/2-1))';
 f=range*fs/N;

 figure(1), clf
 phi = (sqrt(5)+1)/2;
 axes('position', [.125 .15 .75 phi-1])

 ax = [0 100 3e-4 3e2];
 mygridsemycoarse(ax, [0:10:ax(2)], [1e-3 1e-2 1e-1 1e0 1e1 1e2])
 for k=1:n
    plot(freq(k)*[1 1], ax(3:4),'k-','linewidth',1.5)
 end
 semilogy(f,abs(w(range,1)))
 xlabel('f (Hz)')
 axis(ax)
 set(gca,'xtick',[0:10:ax(2)])
 print -depsc2 bead3a