function mfta(rows, cols)
%MFTA(rows, cols)   Multi-Frequency Transform Analysis
%
%  rows (1): how many rows of graphs to plot
%  cols (1): how many columns of graphs to plot
%
%Examples:
%  mfta
%  mfta(5);
%  mfta(5,6);

% PRODUCE MULTIPLE PLOTS FROM RR DATA (script emailed from Thanos to Bob 8/1/00;
% modified a fair bit; later modified by Brandon Beck 5/11/02):
% -------------------------------------------------------------------------------
% RR intervals and instantaneous heart rate graphs

% Tests on the input arguments
if (nargin==0) rows=1; end %if
if (nargin<=1) cols=1; end %if

numgraphs = rows * cols;
for n=1:numgraphs

clear global ecgfullfilename cdfullfilename rrfullfilename hrfullfilename ecgonlyfilename;
try
  [ecgfilename, cdfilename, rrfilename, hrfilename, ecgfilenameonly] = chooser('choosenew');
catch
  break;
end
  
fs = textread(cdfilename, '%f', 1);
[rt rr] = textread(rrfilename, '%f %f', -1, 'headerlines', 1);
[hrt hr] = textread(hrfilename, '%f %f', -1, 'headerlines', 1);
rrsize = size(rr,1);

% ---------------------------------------------------------------------------
% WELCH'S PSD

figure(1); subplot(rows, cols, n);
nfft = 4*1024;
lwdth=1;
fsize=6;
fweight='b';
[p,f] = psd(hr, nfft, fs, kaiser(nfft/4,5), 0 , 0.99,'mean');
[p2,f2] = psd(hr, nfft, fs, hanning(nfft/4), 0 , 0.99,'mean');
plot(f, sqrt(p),'k', 'linewidth', lwdth);
set(gca, 'xlim', [0 10]); set(gca, 'ylim', [0 15]); % z;
hold on
plot(f2, sqrt(p2), 'k--', 'linewidth', lwdth); set(gca, 'fontsize', 8);
set(gca, 'xlim', [0 10]); set(gca, 'ylim', [0 15]); % z;
hold off
xlabel('Frequency (Hz)', 'fontsize', fsize, 'fontweight', fweight );
ylabel('Sqrt(Power)', 'fontsize', fsize, 'fontweight', fweight );
title('Welch modified periodogram', 'fontsize', 8, 'fontweight', 'b')
legend('Kaiser \beta = 4', 'Hanning');
totalpower = (trapz(f,p.^2) + trapz(f2,p2.^2))/2;
text(.5,14,['Total Power = ' num2str(totalpower)], 'fontsize', 8);
text(10.25,13,['Filename = ' ecgfilenameonly], 'fontsize', 8, 'rotation', -90);

% ---------------------------------------------------------------------------
% Yule-Walker AR method
figure(2); subplot(rows, cols, n);
nfft = 1024;  
np = 15; [p,f] = pyulear( hr-median(hr), np, nfft, fs);
np = 15; [p2,f2] = pburg( hr-median(hr), np, nfft, fs);
plot(f, sqrt(p),'k', 'linewidth', lwdth); set(gca, 'fontsize', 8);
set(gca, 'xlim', [0 10]); set(gca, 'ylim', [0 15]); %z;
hold on
plot(f2, sqrt(p2), 'k--', 'linewidth', lwdth); set(gca, 'fontsize', 8);
set(gca, 'xlim', [0 10]); set(gca, 'ylim', [0 15]); %z;
hold off
xlabel('Frequency (Hz)', 'fontsize', fsize, 'fontweight', fweight );
ylabel('Sqrt(Power)', 'fontsize', fsize, 'fontweight', fweight );
title('Autoregressive Methods', 'fontsize', 8, 'fontweight', 'b')
legend('Yule-Walker', 'Burg');
totalpower = (trapz(f,p.^2) + trapz(f2,p2.^2))/2;
text(.5,14,['Total Power = ' num2str(totalpower)], 'fontsize', 8);
text(10.25,13,['Filename = ' ecgfilenameonly], 'fontsize', 8, 'rotation', -90);

% ---------------------------------------------------------------------------
% MUSIC

figure(3); subplot(rows, cols, n);
nfft = 2*1024; [p,f] = pmusic( hr-median(hr), 20, nfft, fs ); 
plot(f, sqrt(p),'k', 'linewidth', lwdth); set(gca, 'fontsize', 8);
set(gca, 'xlim', [0 40]); set(gca, 'ylim', [0 60]); %z;
xlabel('Frequency (Hz)', 'fontsize', fsize, 'fontweight', fweight );
ylabel('Sqrt(Power)', 'fontsize', fsize, 'fontweight', fweight );
title('Multiple Signal Classification Method (MUSIC)', 'fontsize', 8, 'fontweight', 'b')
totalpower = trapz(f,p.^2);
text(13,56,['Total Power = ' num2str(totalpower)], 'fontsize', 8);
text(41.5,59,['Filename = ' ecgfilenameonly], 'fontsize', 8, 'rotation', -90);
%text(3,14,['Total Power = ' num2str(totalpower)], 'fontsize', 8);
%text(10.25,13,['Filename = ' ecgfilenameonly], 'fontsize', 8, 'rotation', -90);

end %for