CEN/EEN330 Course Project Title: Signal Mixing and Source Separation Tools: MATLAB
__MACOSX/._pca_ica
pca_ica/centerRows.m
function [Zc, mu] = centerRows(Z) % % Syntax: [Zc, mu] = centerRows(Z); % % Inputs: Z is an (d x n) matrix containing n samples of a % d-dimensional random vector % % Outputs: Zc is the centered version of Z % % mu is the (d x 1) sample mean of Z % % Description: Returns the centered (zero mean) version of the input data % % Note: Z = Zc + repmat(mu,1,n) % % Author: Brian Moore % brimoor@umich.edu % % Date: November 1, 2016 % % Compute sample mean mu = mean(Z,2); % Subtract mean Zc = bsxfun(@minus,Z,mu);
__MACOSX/pca_ica/._centerRows.m
pca_ica/demo_PCA1.m
% % Name: demo_PCA1.m % % Description: Generates multivariate Gaussian samples and applies PCA to % extract the dimensions of maximum variance and projects the % samples onto them % % Author: Brian Moore % brimoor@umich.edu % % Date: April 26, 2015 % November 12, 2016 % rng(1); % Knobs n = 1000; % Number of samples d = 3; % Sample dimension r = 2; % Number of principal components % Generate Gaussian data MU = 10 * rand(d,1); sigma = (2 * randi([0 1],d) - 1) .* rand(d); SIGMA = 3 * (sigma * sigma'); Z = mvg(MU,SIGMA,n); % Perform PCA [Zpca, ~, ~, eigVecs] = PCA(Z,r); % Plot samples + scaled eigenvectors figure; x = @(i,j) MU(i) + [0 eigVecs(i,j)]; if (d == 2) % Plot 2D data plot(Z(1,:),Z(2,:),'g+'); hold on; for j = 1:r plot(x(1,j),x(2,j),'b','Linewidth',3); end else % Plot first 3 dimensions plot3(Z(1,:),Z(2,:),Z(3,:),'g+'); hold on; for j = 1:r plot3(x(1,j),x(2,j),x(3,j),'b','Linewidth',3); end end title(sprintf('First %i principal directions of %iD Gaussian samples',min(3,r),d)); grid on; set(gca,'DataAspectRatio',[1 1 1]); % Plot principal componenets figure; if (r == 2) % Plot 2D data plot(Zpca(1,:),Zpca(2,:),'r+'); else % Plot first 3 dimensions plot3(Zpca(1,:),Zpca(2,:),Zpca(3,:),'r+'); zlabel('Principal direction #3'); end title(sprintf('Zpca(1:%i,:)',min(3,r))); xlabel('Principal direction #1'); ylabel('Principal direction #2'); grid on; set(gca,'DataAspectRatio',[1 1 1]);
__MACOSX/pca_ica/._demo_PCA1.m
pca_ica/normalizeAudio.m
function Zn = normalizeAudio(Z) % Syntax: Zn = normalizeAudio(Z); % Normalize rows to [-1, 1] Zn = bsxfun(@rdivide,Z,max(abs(Z),[],2));
__MACOSX/pca_ica/._normalizeAudio.m
pca_ica/.DS_Store
__MACOSX/pca_ica/._.DS_Store
pca_ica/mvg.m
function Z = mvg(MU,SIGMA,n) % % Syntax: Z = mvg(MU,SIGMA); % Z = mvg(MU,SIGMA,n); % % Inputs: MU is the desired (d x 1) mean vector % % SIGMA is the desired (d x d) covariance matrix % % [OPTIONAL] n is the desired number of samples. The default % value is n = 1 % % Outputs: Z is a (d x n) matrix containing n multivariate Gaussian % samples % % Description: Generates n samples from a d-dimensional multivariate % Gaussian (MVG) random variable with mean vector MU and % covariance matrix SIGMA % % Notes: (a) mean(Z,2) ~ MU % (b) cov(Z') ~ SIGMA % % Author: Brian Moore % brimoor@umich.edu % % Date: November 1, 2016 % % Parse inputs if ~exist('n','var') || isempty(n) % Default # samples n = 1; end % Generate samples d = size(SIGMA,1); L = chol(SIGMA,'lower'); Z = bsxfun(@plus,MU,L * randn(d,n));
__MACOSX/pca_ica/._mvg.m
pca_ica/demo_PCA2.m
% % Name: demo_PCA2.m % % Description: Demonstrates the performance of PCA by decomposing % correlated multivariate Gaussian samples into % their principal components % % When d == r, the PCA reconstruction is *exact*. Thus % PCA effectively transforms the input data into % orthogonal, maximal variance components while preserving % information % % Author: Brian Moore % brimoor@umich.edu % % Date: April 26, 2015 % November 12, 2016 % rng(42); % Knobs n = 100; % Number of samples d = 3; % Sample dimension r = 2; % Number of principal components % Generate Gaussian data MU = 10 * rand(d,1); sigma = (2 * randi([0 1],d) - 1) .* rand(d); SIGMA = 3 * (sigma * sigma'); Z = mvg(MU,SIGMA,n); % Perform PCA [Zpca, U, mu] = PCA(Z,r); Zr = U * Zpca + repmat(mu,1,n); % Plot principal components figure; for i = 1:r subplot(r,1,i); plot(Zpca(i,:),'b'); grid on; ylabel(sprintf('Zpca(%i,:)',i)); end subplot(r,1,1); title('Principal Components'); % Plot r-dimensional approximations figure; for i = 1:d subplot(d,1,i); hold on; p1 = plot(Z(i,:),'--r'); p2 = plot(Zr(i,:),'-.b'); grid on; ylabel(sprintf('Z(%i,:)',i)); end subplot(d,1,1); title(sprintf('%iD PCA approximation of %iD data',r,d)); legend([p1 p2],'Z','Zr');
__MACOSX/pca_ica/._demo_PCA2.m
pca_ica/fastICA.m
function [Zica, W, T, mu] = fastICA(Z,r,type,flag) % % Syntax: Zica = fastICA(Z,r); % Zica = fastICA(Z,r,type); % Zica = fastICA(Z,r,type,flag); % [Zica, W, T, mu] = fastICA(Z,r); % [Zica, W, T, mu] = fastICA(Z,r,type); % [Zica, W, T, mu] = fastICA(Z,r,type,flag); % % Inputs: Z is an d x n matrix containing n samples of d-dimensional % data % % r is the number of independent components to compute % % [OPTIONAL] type = {'kurtosis','negentropy'} specifies % which flavor of non-Gaussianity to maximize. The default % value is type = 'kurtosis' % % [OPTIONAL] flag determines what status updates to print % to the command window. The choices are % % flag = 0: no printing % flag = 1: print iteration status % % Outputs: Zica is an r x n matrix containing the r independent % components - scaled to variance 1 - of the input samples % % W and T are the ICA transformation matrices such that % Zr = T \ W' * Zica + repmat(mu,1,n); % is the r-dimensional ICA approximation of Z % % mu is the d x 1 sample mean of Z % % Description: Performs independent component analysis (ICA) on the input % data using the Fast ICA algorithm % % Reference: Hyvärinen, Aapo, and Erkki Oja. "Independent component % analysis: algorithms and applications." Neural networks % 13.4 (2000): 411-430 % % Author: Brian Moore % brimoor@umich.edu % % Date: April 26, 2015 % November 12, 2016 % May 4, 2018 % % Constants TOL = 1e-6; % Convergence criteria MAX_ITERS = 100; % Max # iterations % Parse inputs if ~exist('flag','var') || isempty(flag) % Default display flag flag = 1; end if ~exist('type','var') || isempty(type) % Default type type = 'kurtosis'; end n = size(Z,2); % Set algorithm type if strncmpi(type,'kurtosis',1) % Kurtosis USE_KURTOSIS = true; algoStr = 'kurtosis'; elseif strncmpi(type,'negentropy',1) % Negentropy USE_KURTOSIS = false; algoStr = 'negentropy'; else % Unsupported type error('Unsupported type ''%s''',type); end % Center and whiten data [Zc, mu] = centerRows(Z); [Zcw, T] = whitenRows(Zc); % Normalize rows to unit norm normRows = @(X) bsxfun(@rdivide,X,sqrt(sum(X.^2,2))); % Perform Fast ICA if flag % Prepare status updates fmt = sprintf('%%0%dd',ceil(log10(MAX_ITERS + 1))); str = sprintf('Iter %s: max(1 - |<w%s, w%s>|) = %%.4g\\n',fmt,fmt,fmt); fprintf('***** Fast ICA (%s) *****\n',algoStr); end W = normRows(rand(r,size(Z,1))); % Random initial weights k = 0; delta = inf; while delta > TOL && k < MAX_ITERS k = k + 1; % Update weights Wlast = W; % Save last weights Sk = W * Zcw; if USE_KURTOSIS % Kurtosis G = 4 * Sk.^3; Gp = 12 * Sk.^2; else % Negentropy G = Sk .* exp(-0.5 * Sk.^2); Gp = (1 - Sk.^2) .* exp(-0.5 * Sk.^2); end W = (G * Zcw') / n - bsxfun(@times,mean(Gp,2),W); W = normRows(W); % Decorrelate weights [U, S, ~] = svd(W,'econ'); W = U * diag(1 ./ diag(S)) * U' * W; % Update convergence criteria delta = max(1 - abs(dot(W,Wlast,2))); if flag fprintf(str,k,k,k - 1,delta); end end if flag fprintf('\n'); end % Independent components Zica = W * Zcw;
__MACOSX/pca_ica/._fastICA.m
pca_ica/whitenRows.m
function [Zw, T] = whitenRows(Z) % % Syntax: [Zw, T] = whitenRows(Z); % % Inputs: Z is an (d x n) matrix containing n samples of a % d-dimensional random vector % % Outputs: Zw is the whitened version of Z % % T is the (d x d) whitening transformation of Z % % Description: Returns the whitened (identity covariance) version of the % input data % % Notes: (a) Must have n >= d to fully whitenRows Z % % (b) Z = T \ Zcw % % Author: Brian Moore % brimoor@umich.edu % % Date: November 1, 2016 % % Compute sample covariance R = cov(Z'); % Whiten data [U, S, ~] = svd(R,'econ'); T = U * diag(1 ./ sqrt(diag(S))) * U'; Zw = T * Z;
__MACOSX/pca_ica/._whitenRows.m
pca_ica/randomMixingMatrix.m
function A = randomMixingMatrix(d,p) % Syntax: A = randomMixingMatrix(d,p); A = 0.25 + rand(d,p); A = bsxfun(@rdivide,A,sum(A,2));
__MACOSX/pca_ica/._randomMixingMatrix.m
pca_ica/PlayAudio.m
classdef PlayAudio < handle % % Class for playing and visualizing audio files % % Brian Moore % brimoor@umich.edu % % November 10, 2016 % % % Private constants % properties (GetAccess = private, Constant = true) % Audio FPS = 15; % Graphics playback, frames/sec LATENCY_MS = 500; % Audioplayer latency, in ms WINDOW_MS = 50; % Spectrogram window width, in ms % GUI GUI_NAME = 'Audio Player'; % GUI name Y_PAD = 0.05; % Ampitude (y-axis) padding LINEWIDTH = 1; % Status linewidth FONT_SIZE = 10; % Font size AXIS_GAP = 40; % Axis gap, in pixels LIST_GAP = 20; % List gap, in pixels TOP_GAP = 35; % Top gap, in pixels LEFT_GAP = 60; % Left gap, in pixels LIST_WIDTH = 150; % List width, in pixels BUTTON_HEIGHT = 20; % Button height, in pixels % Colors WAVE_COLOR = [93, 147, 191] / 255; % Waveform color STATUS_COLOR = [0, 0, 0] / 255; % Waveform status color PLAY_COLOR = [252, 252, 252] / 255; % Active color STOP_COLOR = [204, 51, 51] / 255; % Stop color end % % Public GetAccess properties % properties (GetAccess = public, SetAccess = private) isPlaying = false; % Playback status currIdx = nan; % Current track index end % % Private properties % properties (Access = public) % Data audio; % Audio data % Graphics fig; % Figure handle ax1; % Waveform axis handle ax2; % Spectrogram axis handle ax3; % Track list title axis handle listh; % Audio list handle playh; % Play button handle ph1; % Waveform slider handle ph2; % Spectrofram slider handle xLim; % Axes x-limits yLim1; % Waveform y-limits yLim2; % Spectrogram y-limits end % % Public methods % methods (Access = public) % % Constructor % function this = PlayAudio(audio) % Syntax: player = PlayAudio(); % player = PlayAudio(audio); % Parse inputs this.audio = audio; % Initialize GUI this.InitializeGUI(); end % % Select track % function Select(this,idx) idx = max(1,min(idx,numel(this.audio))); this.SetTrack(idx); end % % Start playback % function Start(this) if ~this.isPlaying resume(this.audio(this.currIdx).player); end this.UpdateButton(); end % % Stop playback % function Stop(this) if this.isPlaying pause(this.audio(this.currIdx).player); end this.UpdateButton(); end % % Toggle playback % function TogglePlayback(this) if this.isPlaying this.Stop(); else this.Start(); end end % % Reset playback % function Reset(this) stop(this.audio(this.currIdx).player); this.UpdateButton(); this.UpdateStatusLine(); end % % Close player % function Close(this) try this.Reset(); catch %#ok % Empty end try delete(this.fig); catch %#ok delete(gcf); end end end % % Private methods % methods (Access = private) % % Handle list selection % function HandleListSelection(this) value = get(this.listh,'Value'); this.SetTrack(value); end % % Handle keypress % function HandleKeyPress(this,e) % Parse keypress keyChar = e.Character; if isempty(keyChar) return; end % Handle keypress switch double(keyChar) case 13 % Enter % Toggle playback this.TogglePlayback(); case {28, 8, 127} % {Left arrow, backspace, delete} % Reset playback this.Reset(); end end % % Set playing status % function SetPlayStatus(this,isPlaying) this.isPlaying = isPlaying; this.UpdateButton(); this.UpdateStatusLine(); end % % Set track % function SetTrack(this,idx) if idx == this.currIdx return; elseif ~isnan(this.currIdx) this.Reset(); this.audio(this.currIdx).player = []; end % Update track this.currIdx = idx; % Get data y = max(-1,min(this.audio(this.currIdx).y,1)); Fs = this.audio(this.currIdx).Fs; Ny = numel(y); t = (1:Ny) / Fs; % Compute spectrogram window = this.WINDOW_MS * (Fs / 1000); [S, F, T] = spectrogram(y,window,[],[],Fs,'yaxis'); P = 10 * log10(abs(S)); % Update axis limits this.xLim = [t(1), t(Ny)]; this.yLim1 = [-1, 1]; this.yLim2 = [min(F), max(F)]; % Update waveform axLim1 = [this.xLim, this.yLim1 + this.Y_PAD * [-1, 1]]; set(this.ph1(1),'XData',t,'YData',y); axis(this.ax1,axLim1); % Update spectrogram axLim2 = [this.xLim, this.yLim2]; set(this.ph1(2),'XData',T, ... 'YData',F, ... 'ZData',zeros(size(P)), ... 'CData',P); axis(this.ax2,axLim2); % Construct audio player player = audioplayer(y,Fs); player.TimerPeriod = round(1000 / this.FPS) / 1000; player.StartFcn = @(s,e)this.SetPlayStatus(true); player.StopFcn = @(s,e)this.SetPlayStatus(false); player.TimerFcn = @(s,e)this.UpdateStatusLine(); this.audio(this.currIdx).player = player; end % % Update track list % function UpdateTrackList(this) set(this.listh,'String',{this.audio.name}); end % % Update button % function UpdateButton(this) if this.isPlaying set(this.playh,'String','Stop' ,... 'BackgroundColor',this.STOP_COLOR, ... 'ForegroundColor',[1, 1, 1]); else set(this.playh,'String','Play', ... 'BackgroundColor',this.PLAY_COLOR, ... 'ForegroundColor',[0, 0, 0]); end end % % Update status line % function UpdateStatusLine(this) % Get current time if this.isPlaying n = this.audio(this.currIdx).player.CurrentSample; Fs = this.audio(this.currIdx).Fs; else n = -1; Fs = 1; end t = (n / Fs) - (this.LATENCY_MS / 1000); % Update waveform status line x1 = [t, t]; y1 = this.yLim1; set(this.ph2(1),'XData',x1,'YData',y1); % Update spectrogram status line x2 = [t, t]; y2 = this.yLim2; set(this.ph2(2),'XData',x2,'YData',y2); end % % Initialize GUI % function InitializeGUI(this) % Setup figure this.fig = figure( ... 'Name',this.GUI_NAME, ... 'MenuBar','none', ... 'NumberTitle','off', ... 'DockControl','off', ... 'KeyPressFcn',@(s,e)this.HandleKeyPress(e), ... 'ResizeFcn',@(s,e)this.ResizeGUI(), ... 'CloseRequestFcn',@(s,e)this.Close(), ... 'Visible','off'); % File menu filem = uimenu(this.fig,'Label','File'); uimenu(filem,'Label','Close', ... 'Callback',@(s,e)this.Close(), ... 'Accelerator','W'); % Play menu playm = uimenu(this.fig,'Label','Play'); uimenu(playm,'Label','Start', ... 'Callback',@(s,e)this.Start(), ... 'Accelerator','1'); uimenu(playm,'Label','Stop', ... 'Callback',@(s,e)this.Stop(), ... 'Accelerator','2'); uimenu(playm,'Label','Reset', ... 'Callback',@(s,e)this.Reset(), ... 'Accelerator','3', ... 'Separator','on'); % Track list this.ax3 = axes('Units','pixels'); title(this.ax3,'Tracks'); axis(this.ax3,'off'); this.listh = uicontrol( ... 'Parent',this.fig, ... 'Style','List', ... 'Value',1, ... 'Max',1, ... 'Enable','on', ... 'Units','pixels', ... 'HorizontalAlignment','left', ... 'FontSize',this.FONT_SIZE, ... 'Callback',@(s,e)this.HandleListSelection(), ... 'KeyPressFcn',@(s,e)this.HandleKeyPress(e)); % Play button this.playh = uicontrol( ... 'Parent',this.fig, ... 'Style','pushbutton', ... 'String','', ... 'Units','pixels', ... 'FontSize',this.FONT_SIZE, ... 'HorizontalAlignment','center', ... 'Callback',@(s,e)this.TogglePlayback()); % Waveform this.ax1 = axes('Units','pixels'); this.ph1(1) = plot( ... this.ax1,nan,nan,'-','Color',this.WAVE_COLOR); ylabel(this.ax1,'Amplitude','FontSize',this.FONT_SIZE); title(this.ax1,'Waveform','FontSize',this.FONT_SIZE); set(this.ax1,'FontSize',this.FONT_SIZE); axis(this.ax1,'manual'); hold(this.ax1,'on'); % Spectrogram this.ax2 = axes('Units','pixels'); this.ph1(2) = pcolor(this.ax2,nan(2)); xlabel(this.ax2,'Time (s)','FontSize',this.FONT_SIZE); ylabel(this.ax2,'Frequency (Hz)','FontSize',this.FONT_SIZE); title(this.ax2,'Spectrogram','FontSize',this.FONT_SIZE); set(this.ax2,'FontSize',this.FONT_SIZE); shading(this.ax2,'flat'); hold(this.ax2,'on'); % Waveform status line this.ph2(1) = plot( ... nan,nan,'-', ... 'Color',this.STATUS_COLOR, ... 'Linewidth',this.LINEWIDTH, ... 'Parent',this.ax1); % Spectrogram waveform status line this.ph2(2) = plot( ... nan,nan,'-', ... 'Color',this.STATUS_COLOR, ... 'Linewidth',this.LINEWIDTH, ... 'Parent',this.ax2); % Initialize GUI this.UpdateButton(); this.UpdateTrackList(); this.SetTrack(1); this.ResizeGUI(); set(this.fig,'Visible','on'); end % % Resize GUI % function ResizeGUI(this) % Get constants gapa = this.AXIS_GAP; gapb = this.LIST_GAP; gapt = this.TOP_GAP; gapl = this.LEFT_GAP; w2 = this.LIST_WIDTH; bh = this.BUTTON_HEIGHT; % Get figure dimensions pos = get(this.fig,'Position'); xt = pos(3); yt = pos(4); % Compute GUI element positions w1 = xt - w2 - gapl - 2 * gapb; h1 = 0.5 * (yt - 2 * gapa - gapt); h2 = yt - gapa - gapb - gapt - bh; ax1Pos = [gapl, (2 * gapa + h1), w1, h1]; ax2Pos = [gapl, gapa, w1, h1]; listPos = [(gapl + gapb + w1), (gapa + gapb + bh), w2, h2]; playPos = [(gapl + gapb + w1), gapa, w2, bh]; % Update positions set(this.ax1,'Position',ax1Pos); set(this.ax2,'Position',ax2Pos); set(this.ax3,'Position',listPos); set(this.listh,'Position',listPos); set(this.playh,'Position',playPos); end end end
__MACOSX/pca_ica/._PlayAudio.m
pca_ica/visualizeAudio.m
function visualizeAudio(y,name,Fs) % Syntax: visualizeAudio(y,name); % visualizeAudio(y,name,Fs); % Constants WINDOW_MS = 40; % Spectrogram window width, in ms Y_PAD = 0.05; % Ampitude (y-axis) padding % Parse inputs if ~exist('Fs','var') || isempty(Fs) % Default sampling rate Fs = 8000; end Ny = numel(y); t = (1:Ny) / Fs; axLim = [t(1), t(Ny), -(1 + Y_PAD), (1 + Y_PAD)]; % Compute spectrogram window = WINDOW_MS * (Fs / 1000); [S, F, T] = spectrogram(y,window,[],[],Fs,'yaxis'); P = 10 * log10(abs(S)); size(P) size(F) size(T) % Plot waveform subplot(2,1,1); plot(t,y); ax1 = gca(); axis(ax1,axLim); axis(ax1,'manual'); hold(ax1,'on'); xlabel('Time (s)'); ylabel('Amplitude'); title('Waveform'); % Plot spectrogram subplot(2,1,2); pcolor(T,F,P); shading flat; xlabel('Time (s)'); ylabel('Frequency (Hz)'); title('Spectrogram'); % Set figure name fig = gcf(); set(fig,'Name',name); % Play sound sound(y);
__MACOSX/pca_ica/._visualizeAudio.m
pca_ica/loadAudio.m
function Z = loadAudio(paths) % Syntax: Z = loadAudio(paths); % Parse inputs p = numel(paths); % Load audio audio = cell(1,p); for i = 1:p path = paths{i}; [~, ~, ext] = fileparts(path); if strcmpi(ext,'.mat') % Load from .mat file data = load(path); fields = fieldnames(data); audio{i} = data.(fields{1}); else % Load from .wav file try % Old MATLAB versions audio{i} = wavread(path,'double'); catch %#ok % New MATLAB versions audio{i} = audioread(path,'double'); end end end % Determine sample length nSamples = cellfun(@numel,audio); n = min(nSamples); % Construct audio matrix Z = zeros(p,n); for i = 1:p ni = nSamples(i); gapL = floor(0.5 * (ni - n)); gapR = ceil( 0.5 * (ni - n)); Z(i,:) = audio{i}((gapL + 1):(ni - gapR)); end % Normalize audio Z = normalizeAudio(Z);
__MACOSX/pca_ica/._loadAudio.m
pca_ica/kICA.m
function [Zica, W, T, mu] = kICA(Z,r) % % Syntax: Zica = kICA(Z,r); % [Zica, W, T, mu] = kICA(Z,r); % % Inputs: Z is an d x n matrix containing n samples of d-dimensional % data % % r is the number of independent components to compute % % Outputs: Zica is an r x n matrix containing the r independent % components - scaled to variance 1 - of the input samples % % W and T are the ICA transformation matrices such that % Zr = T \ W' * Zica + repmat(mu,1,n); % is the r-dimensional ICA approximation of Z % % mu is the d x 1 sample mean of Z % % Description: Performs independent component analysis (ICA) on the input % data using the max-kurtosis ICA algorithm % % Reference: http://www.cs.nyu.edu/~roweis/kica.html % % Author: Brian Moore % brimoor@umich.edu % % Date: November 12, 2016 % % Center and whiten data mu = mean(Z,2); T = sqrtm(inv(cov(Z'))); Zcw = T * bsxfun(@minus,Z,mu); % Max-kurtosis ICA [W, ~, ~] = svd(bsxfun(@times,sum(Zcw.^2,1),Zcw) * Zcw'); Zica = W(1:r,:) * Zcw;
__MACOSX/pca_ica/._kICA.m
pca_ica/demo_ICA.m
% % Name: demo_ICA.m % % Description: Demonstrates the use of ICA for source separation % % Author: Brian Moore % brimoor@umich.edu % % Date: November 12, 2016 % %% ICA demo (signals) rng(42); % Knobs n = 1000; % # samples T = [3, 4, 5]; % # periods for each signal SNR = 50; % Signal SNR d = 3; % # mixed observations r = 3; % # independent/principal components % Generate ground truth t = @(n,T) linspace(0,1,n) * 2 * pi * T; Ztrue(1,:) = sin(t(n,T(1))); % Sinusoid Ztrue(2,:) = sign(sin(t(n,T(2)))); % Square Ztrue(3,:) = sawtooth(t(n,T(3))); % Sawtooth % Add some noise to make the signals "look" interesting sigma = @(SNR,X) exp(-SNR / 20) * (norm(X(:)) / sqrt(numel(X))); Ztrue = Ztrue + sigma(SNR,Ztrue) * randn(size(Ztrue)); % Generate mixed signals normRows = @(X) bsxfun(@rdivide,X,sum(X,2)); A = normRows(rand(d,3)); Zmixed = A * Ztrue; % Perform Fast ICA Zfica = fastICA(Zmixed,r); % Perform max-kurtosis ICA Zkica = kICA(Zmixed,r); % Perform PCA Zpca = PCA(Zmixed,r); %-------------------------------------------------------------------------- % Plot results %-------------------------------------------------------------------------- cm = hsv(max([3, r, d])); %cm = linspecer(max([3, r, d])); figure(); % Truth subplot(5,1,1); for i = 1:3 plot(Ztrue(i,:),'-','Color',cm(i,:)); hold on; end title('True signals'); axis tight; % Observations subplot(5,1,2); for i = 1:d plot(Zmixed(i,:),'-','Color',cm(i,:)); hold on; end title('Observed mixed signals'); axis tight; % Fast ICA subplot(5,1,3); for i = 1:r plot(Zfica(i,:),'-','Color',cm(i,:)); hold on; end title('Independent components [Fast ICA]'); axis tight; % Max-kurtosis subplot(5,1,4); for i = 1:r plot(Zkica(i,:),'-','Color',cm(i,:)); hold on; end title('Independent components [max-kurtosis]'); axis tight; % PCA subplot(5,1,5); for i = 1:r plot(Zpca(i,:),'-','Color',cm(i,:)); hold on; end title('Principal components'); axis tight; %-------------------------------------------------------------------------- %% ICA demo (audio) % % Audio descriptions % % source1: siren % source2: news (stocks) % source3: foreign 1 % source4: news (food) % source5: music (classical) % source6: foreign 2 % source7: music (opera) % source7: foreign 3 % source9: music (pop) % % voice1: talking (linear algebra) % voice2: talking (sports) % rng(42); % Knobs Fs = 8000; % Sampling rate of input audio paths = {'audio/source2.wav','audio/source3.wav','audio/source5.wav'}; %paths = {'audio/voice1.mat','audio/voice2.mat'}; [p, d, r] = deal(numel(paths)); A = randomMixingMatrix(d,p); %A = [0.6, 0.4; 0.4, 0.6]; % Load audio Ztrue = loadAudio(paths); % Generate mixed signals Zmixed = normalizeAudio(A * Ztrue); % Fast ICA (negentropy) Zica1 = normalizeAudio(fastICA(Zmixed,r,'negentropy')); % Fast ICA (kurtosis) Zica2 = normalizeAudio(fastICA(Zmixed,r,'kurtosis')); % Max-kurtosis ICA Zica3 = normalizeAudio(kICA(Zmixed,r)); %-------------------------------------------------------------------------- % Plot results %-------------------------------------------------------------------------- audio = []; % True waveforms for i = 1:p audio(end + 1).y = Ztrue(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('true - %d',i); end % Mixed waveforms for i = 1:d audio(end + 1).y = Zmixed(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('mixed - %d',i); end % Fast ICA (negentropy) waveforms for i = 1:r audio(end + 1).y = Zica1(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('fastICA - neg - %d',i); end % Fast ICA (kurtosis) waveforms for i = 1:r audio(end + 1).y = Zica2(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('fastICA - kur - %d',i); end % Max-kurtosis ICA waveforms for i = 1:r audio(end + 1).y = Zica3(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('kICA - %d',i); end % Play audio PlayAudio(audio); %-------------------------------------------------------------------------- %% ICA demo (pre-mixed audio) % TODO(brimoor): how to get this to work? % % Audio descriptions % % rsm2_mA/B: music + counting english % rss_mA/B: counting english + counting spanish % rssd_A/B: talking english + talking spanish % rng(42); % Knobs Fs = 16000; % Sampling rate of input audio paths = {'audio/rsm2_mA.wav','audio/rsm2_mB.wav'}; %paths = {'audio/rss_mA.wav','audio/rss_mB.wav'}; %paths = {'audio/rssd_A.wav','audio/rssd_B.wav'}; [d, r] = deal(numel(paths)); r = r + 1; % Load audio Zmixed = loadAudio(paths); % Fast ICA (negentropy) Zica1 = normalizeAudio(fastICA(Zmixed,r,'negentropy')); % Fast ICA (kurtosis) Zica2 = normalizeAudio(fastICA(Zmixed,r,'kurtosis')); % Max-kurtosis ICA Zica3 = normalizeAudio(kICA(Zmixed,min(r,d))); %-------------------------------------------------------------------------- % Plot results %-------------------------------------------------------------------------- audio = []; % Mixed waveforms for i = 1:d audio(end + 1).y = Zmixed(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('mixed - %d',i); end % Fast ICA (negentropy) waveforms for i = 1:r audio(end + 1).y = Zica1(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('fastICA - neg - %d',i); end % Fast ICA (kurtosis) waveforms for i = 1:r audio(end + 1).y = Zica2(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('fastICA - kur - %d',i); end % Max-kurtosis ICA waveforms for i = 1:min(r,d) audio(end + 1).y = Zica3(i,:); %#ok audio(end).Fs = Fs; audio(end).name = sprintf('kICA - %d',i); end % Play audio PlayAudio(audio); %--------------------------------------------------------------------------
__MACOSX/pca_ica/._demo_ICA.m
pca_ica/PCA.m
function [Zpca, U, mu, eigVecs] = PCA(Z,r) % % Syntax: Zpca = PCA(Z,r); % [Zpca, U, mu] = PCA(Z,r); % [Zpca, U, mu, eigVecs] = PCA(Z,r); % % Inputs: Z is an d x n matrix containing n samples of d-dimensional % data % % r is the number of principal components to compute % % Outputs: Zpca is an r x n matrix containing the r principal % components - scaled to variance 1 - of the input samples % % U is a d x r matrix of coefficients such that % Zr = U * Zpca + repmat(mu,1,n); % is the r-dimensional PCA approximation of Z % % mu is the d x 1 sample mean of Z % % eigVecs is a d x r matrix containing the scaled % eigenvectors of the sample covariance of Z % % Description: Performs principal component analysis (PCA) on the input % data % % Author: Brian Moore % brimoor@umich.edu % % Date: April 26, 2015 % November 7, 2016 % % Center data [Zc, mu] = centerRows(Z); % Compute truncated SVD %[U, S, V] = svds(Zc,r); % Equivalent, but usually slower than svd() [U, S, V] = svd(Zc,'econ'); U = U(:,1:r); S = S(1:r,1:r); V = V(:,1:r); % Compute principal components Zpca = S * V'; %Zpca = U' * Zc; % Equivalent but slower if nargout >= 4 % Scaled eigenvectors eigVecs = bsxfun(@times,U,diag(S)' / sqrt(size(Z,2))); end