% Experiments conducted 16 Feb 2006 by Sean Hardesty and Erin Noel
% One bead, pluck without lights

 n = 1;

 m = 0.125;            % .5 oz
 sigma = 3.04*1.5;     % 5.5 lb/mv*mv
 
% units in cm -- reported by Scott
 ell0 = 24.5;   % cm
 ell1 = 55;     % cm

 a1 = sigma/2*(1/ell0+1/ell1);
 

 C = .5*m*eye(1);
 A = [ a1 ];

 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 = 'record002.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 140 1e-0 1e4];
 mygridsemycoarse(ax, [0:20:ax(2)], [1e1 1e2 1e3])
 for k=1:n
    plot(freq(k)*[1 1], ax(3:4),'k--','linewidth',1.5)
 end
 semilogy(f,abs(w(range,1)),'linewidth',1.5)
 xlabel('f (Hz)')
 title(sprintf('%s',filename))
 axis(ax)
 set(gca,'xtick',[0:20:ax(2)])
 print -depsc2 bead1b

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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