%以下為主程序，主要為原始信號的產(chǎn)生，觀察信號和解混信號的作圖
clear all;clc;
N=200;n=1:N;%N為采樣點(diǎn)數(shù)
s1=2*sin(0.02*pi*n);%正弦信號
t=1:N;s2=2*square(100*t,50);%方波信號
a=linspace(1,-1,25);s3=2*[a,a,a,a,a,a,a,a];%鋸齒信號
s4=rand(1,N);%隨機(jī)噪聲
S=[s1;s2;s3;s4];%信號組成4*N
A=rand(4,4);
X=A*S;%觀察信號
%源信號波形圖
figure(1);subplot(4,1,1);plot(s1);axis([0N -5,5]);title('源信號');
subplot(4,1,2);plot(s2);axis([0N -5,5]);
subplot(4,1,3);plot(s3);axis([0N -5,5]);
subplot(4,1,4);plot(s4);xlabel('Time/ms');
%觀察信號（混合信號）波形圖
figure(2);subplot(4,1,1);plot(X(1,:));title('觀察信號（混合信號）');
subplot(4,1,2);plot(X(2,:));
subplot(4,1,3);plot(X(3,:));subplot(4,1,4);plot(X(4,:));
Z=ICA(X);
figure(3);subplot(4,1,1);plot(Z(1,:));title('解混后的信號');
subplot(4,1,2);plot(Z(2,:));
subplot(4,1,3);plot(Z(3,:));
subplot(4,1,4);plot(Z(4,:));xlabel('Time/ms');
function Z=ICA(X)
%-----------去均值---------
[M,T] = size(X); %獲取輸入矩陣的行/列數(shù)，行數(shù)為觀測數(shù)據(jù)的數(shù)目钧舌，列數(shù)為采樣點(diǎn)數(shù)
average=mean(X')'; %均值
for i=1:M
X(i,:)=X(i,:)-average(i)*ones(1,T);
end
%---------白化/球化------
Cx =cov(X',1); %計(jì)算協(xié)方差矩陣Cx
[eigvector,eigvalue]= eig(Cx); %計(jì)算Cx的特征值和特征向量
W=eigvalue^(-1/2)*eigvector'; %白化矩陣
Z=W*X; %正交矩陣
%----------迭代-------
Maxcount=10000; %最大迭代次數(shù)
Critical=0.00001; %判斷是否收斂
m=M; %需要估計(jì)的分量的個(gè)數(shù)
W=rand(m);
for n=1:m
WP=W(:,n); %初始權(quán)矢量（任意）
% Y=WP'*Z;
% G=Y.^3;%G為非線性函數(shù)，可取y^3等
% GG=3*Y.^2; %G的導(dǎo)數(shù)
count=0;
LastWP=zeros(m,1);
W(:,n)=W(:,n)/norm(W(:,n));
while (abs(WP-LastWP)& abs(WP+LastWP)>Critical)
count=count+1; %迭代次數(shù)
LastWP=WP; %上次迭代的值
% WP=1/T*Z*((LastWP'*Z).^3)'-3*LastWP;
for i=1:m
WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i);
end
WPP=zeros(m,1);
for j=1:n-1
WPP=WPP+(WP'*W(:,j))*W(:,j);
end
WP=WP-WPP;
WP=WP/(norm(WP));
if count==Maxcount
fprintf('未找到相應(yīng)的信號');
return;
end
end
W(:,n)=WP;
end
Z=W'*Z;
% ---------------------
% 作者：whiteinblue
% 來源：CSDN
% 原文：https://blog.csdn.net/whiteinblue/article/details/36366817
% 版權(quán)聲明：本文為博主原創(chuàng)文章络拌，轉(zhuǎn)載請附上博文鏈接乃坤！
