Commit da186c7c authored by tllado's avatar tllado

Added HW2 code

parent 3bae513c
code @ 365ba3cd
Subproject commit 365ba3cd2caa6ee07dd03f0c10ecbc7d1a257182
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hw2ICA
% Mix and then separate audio signals
% written by Travis Llado
% last edited 2017.09.29
% travisllado@utexas.edu
% This program was written using MATLAB v2017a
% Backwards compatibility is not guaranteed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear variables
dataFilename = 'sounds.mat'; % Given, do not change
samplingRate = 11250; % Given, do not change
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% User Settings
% Experiment Parameters
sourcesToUse = [1 4 5 2 3]; % source signals to mix, order 14523 is preferred
numMics = length(sourcesToUse); % # mixed signals to analyze
roomSize = 10; % size of "space" over which signals and microphones are spread, aritrary units
% Solver Parameters
sampleStart = 10001; % # of first data point to use
sampleLength = 500; % # of data points to use
eta = 1e-3; % learning rate
maxIterations = 1e7; % max times to run solver
endCondition = 1e-10; % limit for max(abs(dW))
readFrequency = 1e4; % frequency at which we record error values for analysis
% Options
drawRoom = 1; % Draw map of the source/microphone layout
drawSignals = 0; % Draw plots of the source and mic signals
drawConvergence = 1; % Draw plot of the convergence process
drawSmallComparison = 0; % Plot sources and reconstructions only
drawBigComparison = 1; % Plot sources, microphones, and reconstructions
playMixes = 0; % Play original and reconstructed signals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assemble test data
numSources = length(sourcesToUse);
% Check for solvability
if numSources > numMics
warning('It is strongly recommended that you use more microphones than sources');
response = input(['If you want to continue down this ill advised path, type ''I do what I want''' newline], 's');
if ~strcmp(response, 'I do what I want')
return
end
end
% Retrieve test data
data = open(dataFilename);
sourceSignals = data.sounds(sourcesToUse, sampleStart:sampleStart + sampleLength - 1);
numSamples = length(sourceSignals);
% Place sources within room
sourceCoordinates = zeros(2, numSources);
for sourceNum = 1:numSources
sourceCoordinates(:, sourceNum) = (rand(1, 2) - 0.5)*roomSize;
end
% Place microphones within room
micCoordinates = zeros(2, numMics);
for micNum = 1:numMics
micCoordinates(:, micNum) = (rand(1, 2) - 0.5)*roomSize;
end
% Plot room layout
if drawRoom
figure
plot(sourceCoordinates(1, :), sourceCoordinates(2, :), '+', micCoordinates(1, :), micCoordinates(2, :), 'o', 'Linewidth', 2)
axis([-roomSize/2 roomSize/2 -roomSize/2 roomSize/2])
legend('source', 'microphone')
for sourceNum = 1:numSources
text(sourceCoordinates(1, sourceNum) + roomSize*0.0125, sourceCoordinates(2, sourceNum), num2str(sourceNum))
end
for micNum = 1:numMics
text(micCoordinates(1, micNum) + roomSize*0.0125, micCoordinates(2, micNum), num2str(micNum))
end
end
% Mix signals received by microphones
micSignals = zeros(numMics, numSamples);
for thisMic = 1:numMics
micDistances = hypot(sourceCoordinates(1, :) - micCoordinates(1, thisMic), sourceCoordinates(2, :) - micCoordinates(2, thisMic));
sourceVolumes = 1./(micDistances.^2);
micSignals(thisMic, :) = ((sourceSignals)'*(sourceVolumes)')';
end
% Scale microphone signals
micSignals = ((micSignals)'./max(abs(micSignals)'))';
% Plot signals for comparison
if drawSignals
figure
for sourceNum = 1:numSources
subplot(numSources, 1, sourceNum)
plot(sourceSignals(sourceNum, :)', 'Linewidth', 2)
ylabel(['Source ' num2str(sourceNum)])
end
figure
for micNum = 1:numMics
subplot(numMics, 1, micNum)
hold on
plot(0)
plot(micSignals(micNum, :)', 'Linewidth', 2)
hold off
ylabel(['Mic ', num2str(micNum)])
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Separate Signals
X = micSignals;
W = rand(numSources, numMics)*2/numSources;
Y = W*X;
numIterations = 0;
currentError = endCondition + 1;
allError = zeros(maxIterations/readFrequency, 1);
timeSpent = mod(now, 1)*24*3600; % start counting time
while (currentError > endCondition) && (numIterations < maxIterations)
numIterations = numIterations + 1;
Y = W*X;
Z = 1./(1 + exp(-Y));
dW = eta*(eye(numSources) + (1 - 2*Z)*(Y)')*W;
W = W + dW;
currentError = norm(dW);
if mod(numIterations, round(readFrequency)) == 0
allError(numIterations/readFrequency) = currentError;
disp(['num iterations = ' num2str(numIterations) ', convergence rate = ' num2str(currentError/endCondition)]);
end
end
timeSpent = mod(now, 1)*24*3600 - timeSpent;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display Results
% Sort signals
matches = zeros(numSources, 1);
used = zeros(numSources, 1);
for inputNum = 1:numSources
minDiff = -1;
for outputNum = 1:numSources
if ~used(outputNum)
crntDiff = mean(abs(sourceSignals(inputNum, :) - Y(outputNum, :)*mean(abs(sourceSignals(inputNum, :)))/mean(abs(Y(outputNum, :)))));
if (minDiff == -1) || (crntDiff < minDiff)
minDiff = crntDiff;
matches(inputNum) = outputNum;
end
end
end
used(matches(inputNum)) = 1;
end
% Scale separated signals
separatedSignals = Y;
for signalNum = 1:numSources
separatedSignals(signalNum, :) = (Y(matches(signalNum), :)/max(abs(Y(matches(signalNum), :))));
sourceSignals(signalNum, :) = (sourceSignals(signalNum, :)/max(abs(sourceSignals(signalNum, :))));
end
% Draw convergence plot
if drawConvergence
figure
plot([1:length(allError)]*readFrequency, log(allError), 'Linewidth', 2)
ylabel('log(norm(dW))')
xlabel('# Iterations')
hold off
end
% Draw originals vs reconstructions
if drawSmallComparison
figure
for signalNum = 1:numSources
subplot(numSources, 2, 2*(signalNum - 1)+1)
plot(sourceSignals(signalNum, :))
ylabel(['Signal ' num2str(signalNum)])
xlabel('Original')
subplot(numSources, 2, 2*signalNum)
hold on
plot(0)
plot(separatedSignals(signalNum, :))
hold off
xlabel('Reconstruction')
end
end
% Draw a ridiculously large plot that shows everything
if drawBigComparison
ySize = max([numSources numMics]);
figure
for n = 1:numSources
subplot(ySize, 3, 3*n - 2)
plot(sourceSignals(n, :), 'Linewidth', 2)
ylabel(['Signal ' num2str(n)])
xlabel('Original')
subplot(ySize, 3, 3*n)
hold on;plot(0);
plot(separatedSignals(n, :), 'Linewidth', 2)
hold off
xlabel('Reconstruction')
end
for n = 1:numMics
subplot(ySize, 3, 3*n - 1)
hold on;plot(0);plot(0);
plot(micSignals(n, :), 'Linewidth', 2)
hold off
ylabel(['Signal ' num2str(n)])
xlabel('Microphone')
end
end
% Play signals
if playMixes
for signalNum = 1:numSources
disp(['Signal ' num2str(signalNum)]);
disp('original')
sound(sourceSignals(signalNum, :), samplingRate);
pause(numSamples/samplingRate)
disp('reconstruction')
sound(separatedSignals(signalNum, :), samplingRate);
pause(numSamples/samplingRate + 0.5)
end
end
disp(['Finished ' num2str(numIterations) ' iterations in ' num2str(timeSpent) ' seconds'])
\ No newline at end of file
% hw2ICA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mix and then separate audio signals
% written by Travis Llado
% last edited 2017.09.29
% travisllado@utexas.edu
% This program was written for MATLAB v2017a
% Program Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear variables
dataFilename = 'sounds.mat'; % Given, do not change
samplingRate = 11250; % Given, do not change
sourcesToUse = [1 5 4 2 3]; % source signals to mix
numSignals = length(sourcesToUse); % # mixed signals to analyze
roomSize = 10; % size of "room" for sources, microphones
sampleStart = 10001; % # of first data point to use
sampleLength = 500; % # of data points to use
eta = 1e-3; % learning rate
numberIterations = 1e7; % max times to run solver
% Create Microphone Signals %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data = open(dataFilename);
sourceSignals = data.sounds(sourcesToUse, sampleStart:sampleStart + sampleLength - 1);
micSignals = zeros(numSignals, sampleLength);
sourceCoordinates = zeros(2, numSignals);
micCoordinates = sourceCoordinates;
for signalNum = 1:numSignals
sourceCoordinates(:, signalNum) = (rand(1, 2) - 0.5)*roomSize;
micCoordinates(:, signalNum) = (rand(1, 2) - 0.5)*roomSize;
end
for thisMic = 1:numSignals
micDistances = hypot(sourceCoordinates(1, :) - micCoordinates(1, thisMic), sourceCoordinates(2, :) - micCoordinates(2, thisMic));
sourceVolumes = 1./(micDistances.^2);
micSignals(thisMic, :) = ((sourceSignals)'*(sourceVolumes)')';
end
% Separate Signals %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = ((micSignals)'./max(abs(micSignals)'))';
W = rand(numSignals, numSignals)*2/numSignals;
Y = W*X;
for iterationNumber = 1:numberIterations
Y = W*X;
Z = 1./(1 + exp(-Y));
dW = eta*(eye(numSignals) + (1 - 2*Z)*(Y)')*W;
W = W + dW;
if mod(iterationNumber, numberIterations/1000) == 0
disp([num2str(iterationNumber*100/numberIterations) '%']);
end
end
% Sort signals %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
matches = zeros(numSignals, 1);
used = eye(numSignals);
for signalNum = 1:numSignals
[~, matches(signalNum)] = min(abs(mean((abs(used*sourceSignals) - abs(Y(signalNum,:)))')));
used(matches(signalNum), matches(signalNum)) = NaN;
end
% Plot results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
for signalNum = 1:numSignals
subplot(numSignals, 3, 3*signalNum - 2)
plot((sourceSignals(signalNum, :)/max(abs(sourceSignals(signalNum, :)))), 'r', 'Linewidth', 2)
ylabel(['Signal ' num2str(signalNum)])
xlabel('Original')
subplot(numSignals, 3, 3*signalNum - 1)
plot(X(signalNum, :), 'g', 'Linewidth', 2)
ylabel(['Signal ' num2str(signalNum)])
xlabel('Microphone')
subplot(numSignals, 3, 3*signalNum)
plot((Y(matches(signalNum), :)/max(abs(Y(matches(signalNum), :)))), 'b', 'Linewidth', 2)
ylabel(['Signal ' num2str(signalNum)])
xlabel('Reconstruction')
end
\ No newline at end of file
D = open('sounds.mat'); % Open source data file
S = D.sounds(1:5, 10001:10500); % Extract source signal data
X = rand(5)*S; % Randomly combine source signals to create microphone signals
W = rand(5, 5)*2/5; % Select random values for initial W
for n = 1:1e7 % Iterate toward signal separation
Y = W*X; % Create current reconstruction
Z = 1./(1 + exp(-Y)); % Compare current reconstruction to expected distribution
dW = 1e-3*(eye(5) + (1 - 2*Z)*(Y)')*W; % Calculate gradient toward correct reconstruction
W = W + dW; % Move toward correct reconstruction
end
for n = 1:5 % Plot originals, reconstructions, and microphones
subplot(5, 3, 3*n - 2); plot(S(n, :), 'r'); xlabel('Original')
subplot(5, 3, 3*n - 1); plot(Y(n, :), 'b'); xlabel('Reconstruction')
subplot(5, 3, 3*n); plot(X(n, :), 'g'); xlabel('Microphone')
end
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment