function [train_coefs,test_coefs,dpb,dpa] = discrpur(n_coef,train_data,train_groups,test_data,test_groups,qmf)
% Use:
%	   [train_coefs,test_coefs,dpb,dpa] = discrpur(n_coef,train_data,train_groups,test_data,test_groups,qmf)
%
% Inputs:   n_coef - how many coefficients to extract
%           train_data, test_data - training and testing data with cols as
%             samples (bins) and rows as observations (trials)
%           train_groups, test_groups - one per observation and starting
%             from 1
%           qmf - quadrature mirror filter to use
% Outputs:
%           train_coefs, test_coefs -  coefficients for each discriminant
%             pursuit basis function
%           dpb - the discrminant pursuit basis functions
%           dpa - the amplitudes for each function
%
% Original code by Jon Buckheit, Dept of Statistics, Stanford University.
% Modified by Mark Laubach, The John B Pierce Laboratory
%
% Original code from JB obtained in 1995 and modified and distributed with
% the JB's permission by ML.
%
% Requires the Wavelab package from Stanford's Dept of Statistics and
% Matlab 5 or greater.
%
% Mark Laubach, September 1, 2003

[n_train, n] = size(train_data);
[n_test, xx] = size(test_data);
groupids = unique(train_groups);
g = length(groupids);
D = log2(n);
n_contr = g*(g-1)/2;

avgscores = zeros(n,g);
pkt_tables = zeros(n_contr,n*(D+1));

for i = 1:g,
	x = train_data(train_groups == groupids(i),:);
   avgscores(:,i) = mean(x)';
end

scale = diag(ones(1,n));

k = 0;
md = zeros(1,n_contr);
for i = 1:(g-1),
	for j = (i+1):g,
		contr = avgscores(:,i)-avgscores(:,j);
		k = k+1;
		pkt = Wpanalysis(scale*contr,D,qmf);
		pkt_tables(k,:) = reshape(pkt,1,n*(D+1));
	end
end

coefs = zeros(1,n_coef);
amp = zeros(1,n_contr);
ind = zeros(1,n_contr);

for i=1:n_coef,

	for j=1:n_contr,
		[amp(j),ind(j)] = max(abs(pkt_tables(j,:)));
	end
	
	[a,l] = max(amp);
	dpa(i) = amp(l);
	coefs(i) = ind(l);
	[d,b,k] =ix2pkt(ind(l),D,n);

	for j=1:n_contr,
		pkt = reshape(pkt_tables(j,:),n,D+1);
		a = pkt(pkt2ix(d,b,k,D,n),d+1);
		dpkt = Wpimpulse(pkt,d,b,k,qmf);
		pkt = pkt - a*dpkt;
		pkt_tables(j,:) = reshape(pkt,1,n*(D+1));
	end
	
end

train_coefs = zeros(n_train,n_coef);
for i=1:n_train,
	signal = train_data(i,:);
	pkt = Wpanalysis(signal,D,qmf);
	train_coefs(i,:) = pkt(coefs);
end

test_coefs = zeros(n_test,n_coef);
for i=1:n_test,
	signal = test_data(i,:);
	pkt = Wpanalysis(signal,D,qmf);
	test_coefs(i,:) = pkt(coefs);
end

for i=1:n_coef
    [d,b,k]= ix2pkt(coefs(i),D,n);	% 
    dpb(i,:) = Makewp(d,b,k,qmf,n);
end