ica_pdf.m

这是盲信号的代码 都已经通过编译了 做这方面的同仁可以参考一下 我觉得蛮惯用的

function ica_pdf()
% mixes sound files
clear all;
M = 2;
s1 = wavread('man.wav');
s2 = wavread('women.wav');
s1=s1';
s2=s2';
maxlen = max( length(s2), length(s1) );
ss(1,:) = [ s1 zeros(1, maxlen - length(s1) ) ];
ss(2,:) = [ s2 zeros(1, maxlen - length(s2) ) ];
s(1,:) = (ss(1,:) - mean( ss(1,:) ) ) ./ max(ss(1,:));
s(2,:) = (ss(2,:) - mean( ss(2,:) ) ) ./ max(ss(2,:));
clear s1 s2;
a11 = .8; a12 = .2;
a21 = .3; a22 = .7;
A = [ [ a11 a12 ]; [ a21 a22 ] ];
x = A * s;
% plots
subplot(2,2,1);
plot(s(1,:));
title('original 1');
subplot(2,2,2);
plot(s(2,:));
title('original 2');
subplot(2,2,3);
plot(x(1,:));
title('mixed 1');
subplot(2,2,4);
plot(x(2,:));
title('mixed 2');
% sounds

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Demixing
%
% for 2 sources
%
% x contains the mixed data in the rows
%
% signal mixture has to be whitenend first in order
% to get good results
%[e,d]=pcamat(x);
%[nv wm dwm ]=whitenv(x,e,d);
N=length( x(1,:) ); % signal length
I = eye(size(x,1));
% take out mean and normalize them to the max value
u(1,:)=( x(1,:) - mean(x(1,:)) ) / max( x(1,:) );
u(2,:)=( x(2,:) - mean(x(2,:)) ) / max( x(2,:) );
W = randn( size(u,1) ); % random start
W_1 = ones( size(u,1) );
% main loop
res = 0;
while (res==0)
% for Super-Gaussian
W = [ I - 2*tanh(u)*u' - u*u' ] * W / N;
W = W / norm(W);
% if vectors of W point in the same direction
% from one iteration to the other, we are done.
d1 = W(1,:)*W_1(1,:)' / ( norm(W(1,:))*norm(W_1(1,:)) );
d2 = W(2,:)*W_1(2,:)' / ( norm(W(2,:))*norm(W_1(2,:)) );
[ d1 d2 ]
res=((d1*d2) > 0.9999); % abort condition
W_1=W;
end;
u = W * x; % demixing
% plots
figure(2);
subplot(2,1,1);
plot(u(1,:));
title('separated 1');
subplot(2,1,2);
plot(u(2,:));
title('separated 2');

% calculate mixing matrix estimate
WAVWRITE(u(1,:),16000,16,'y1.wav');
WAVWRITE(u(2,:),16000,16,'y2.wav');
B=W^(-1);
B/max(max(B))*0.8
% sounds, normalized
u(1,:) = u(1,:) / max(u(1,:) );
u(2,:) = u(2,:) / max(u(2,:) );

%sound(u(1,:),16000)
%disp('press key for next sound');
%pause;
%sound(u(2,:),16000);