function [y, K] = ppursuit(h, eta, tol, mxi, x)
%
% [y, K] = ppursuit(h, eta, tol, mxi, x)
%
% Projection Pursuit (probing to estimate gradient).
% This function uses projection pursuit to demix m
% signal mixtures into m estimated source signals.
%
% h Step size for probing the Kurtosis (K) of the
% demixed signals, relative to the magnitude
% of the current demixing vector. IOW, K is
% probed by looking at K(w + h) in the m
% dimensions of w, one at a time. It should be
% noted that w is a unit vector, so there is
% no need to scale the step size by the current
% size of w.
% eta Distance that demixing vector w is adjusted
% in the direction of the estimated gradient of
% K. IOW, wnext = w + eta*g where g is the
% estimated gradient of K, normalized to 1.
% tol Stopping criterion for gradient ascent. The
% ascent terminates when the relative change in
% K is < tol. (Change in K divided by K).
% mxi Maximum iterations to execute on the ascent
% of K for each recovered signal.
% x (m x n) matrix of signal mixtures. Each row
% is a 1 x n mixture of m source signals. The
% source signals may be anything, and there is
% no assumed relationship between the n samples
% of a given mixture. They may be temporal or
% spatial samples.
% y The (m x n) source estimates
% K The (m x ?) history of the ascent of Kurtosis
% for each recovered source. Having this history
% allows the caller to tune the search parameters.
%
% This version written by Andrew Sackville-West
% for Dr. Schimpf's CSCD309 course
% define our kurtosis function. Because we whiten our data, and set
% the variance to 1, we can use a simplified kurtosis function,
% focusing on the denominator.
k = @(data) mean(data .^ 4) - 3;
% first, we project x onto it's principal components.
% using SVD, we can get a ortho-normal pca basis using just V'
% this allows for easier extraction of the data later
[U S V] = svd(x,0);
% our pca basis, the sound file expressed in an orthonormal basis
z = V';
% determine the number of signals
[signals temp] = size(x);
% loop for each signal
for signal = 1:signals
% set each row to a variance of 1
for row = 1:signals
z(row,:)=z(row,:)./std(z(row,:));
end
% set our starting demixing vector
w = [1,0];
% compute a trial demix
curr_y = w*z;
% setup a matrix of h values for use calculating test w's
hs = diag(ones(1,signals) .* h);
% get the first kurtosis value and set up history for kurtosis values
currK = k(curr_y);
iter = 1;
K(signal,iter) = currK;
deltaK = tol * currK + 1; % force the loop to signal the first time
% loop, incrementing our w vector until we get
% tight enough tolerance, or run out of iterations
% this loop does the probing and pursues the gradient "uphill"
while (abs(deltaK/currK) > tol & iter <= mxi)
% expand w into the right number of rows
test_ws = repmat(w,signals,1);
% add the hs to get w adjusted in each dimension
test_ws = test_ws + hs;
% compute all the test_ys at once
test_ys = test_ws * z;
% calculate the deltaK for each dimension, the gradient
for dim = 1:signals
% calculate the deltaK for the current row of test_y
gradient(dim)= (k(test_ys(dim,:)) - currK)/h;
end % calculating gradient
% modify the gradient by eta
gradient = gradient * eta;
% add the gradient to w and normalize it
w = w + gradient;
w = w./norm(w);
% setup deltaK
deltaK = currK;
% set the new trial kurtosis
currK = k(w*z);
% compute new deltaK
deltaK = deltaK - currK;
% increment our iteration counter
iter += 1;
% store the historic kurtosis values
K(signal,iter) = currK;
end % the pursuit while loop
% calculate the current unmixed signal by projecting z onto the w
% vector
unmix = w * z;
% build the y result
y(signal,:) = unmix;
% calculate a new z for further projection
if(signal ~= signals)
z = z - (w' * unmix);
end
end % main for loop