function analyze_file(filename, f0, Fs)
% function analyze_file(filename, f0, Fs)
% e.g., analyze_file('20hz_60mv.txt', 20)
% f0 = drive frequency

% assign drive frequency variable
f0 = strassign(f0);

% default sampling rate if necessary
try
  Fs;
catch
  Fs = 1000;
  fprintf('assuming data is sampled at 1000 Hz\n');
end

% load and parse data file
g = load(filename);
q = g(:,1);
h = g(:,2);

% Fourier transform data
[hf f] = fft_eht(h, Fs);
[qf f] = fft_eht(q, Fs);

% zoom in on interesting frequencies
cut = f>0.25 & f<65;
%cut = f>f0-1 & f<f0+1;
hf = hf(cut);
qf = qf(cut);
f = f(cut);

% calculate information about peak
zoom = f>=f0-1 & f<=f0+1;
fzoom = f(zoom);
hfzoom = hf(zoom);
qfzoom = qf(zoom);
% top
[qmax, maxidx] =  max(abs(qf(zoom)));
fmaxq = fzoom(maxidx);
fprintf('top:\n');
fprintf('  peak location = %1.2f\n', fmaxq);
fprintf('  peak value = %1.2f\n', qmax);
% bottom
[hmax, maxidx] =  max(abs(hf(zoom)));
fmaxh = fzoom(maxidx);
fprintf('bottom:\n');
fprintf('  peak location = %1.2f\n', fmaxh);
fprintf('  peak value = %1.2f\n', hmax);
% calculate transfer function
T = hmax/qmax;
fprintf('T(f=%1.2f) = %1.2e\n', fmaxh, T);
% calculate FWHM
%fpeakh = fzoom(abs(hf(zoom)) >= hmax/2);
%fpeakq = fzoom(abs(qf(zoom)) >= qmax/2);
hpeakidx = abs(hf(zoom)) >= hmax/2;
qpeakidx = abs(qf(zoom)) >= qmax/2;
T2 = sum(abs(hfzoom(hpeakidx))) / sum(abs(qfzoom(qpeakidx)));
fprintf('T2(f=%1.2f) = %1.2e\n', fmaxh, T2);

% plot results
figure;
subplot(2,1,1);
loglog(f, abs(qf));
hold on;
loglog([f0 f0], [min(abs(qf)) max(abs(qf))], 'r:');
grid on;
axis([f(1) f(end), min(abs(qf)) max(abs(qf))]);
ylabel('top');
pretty;
%
subplot(2,1,2);
loglog(f, abs(hf));
hold on;
loglog([f0 f0], [min(abs(hf)) max(abs(hf))], 'r:');
grid on;
axis([f(1) f(end), min(abs(hf)) max(abs(hf))]);
xlabel('f (Hz)');
ylabel('bottom');
pretty;
print('-dpng', 'spectrum');

return
