function?fushangST1（）
%%?清空環境
clear;clc;
%%?參數設置
fs=1000;
N=100;?%?時間橫坐標軸
t=（0:N-1）/fs;
f1=60;
s1=sin（2*pi*f1*t）;size（s1）?%?正弦信號
s2=2*square（100*t50）;?size（s2）?%方波信號
s3=randn（size（t））;size（s3）?%?隨機信號
%%
%%?主程序
%?輸出原始源信號
subplot（331）plot（s1）title（‘s1‘）axis（[0100-44]）;?%?（331）表示整幅圖分為3行3列共9幅小圖，1表示序號
subplot（332）plot（s2）title（‘s2‘）axis（[0100-44]）;
subplot（333）plot（s3）title（‘s3‘）axis（[0100-44]）;

%?將其組成矩陣
S=[s1;s2;s3];?size（S）;????%?信號個數即為變量數，因此S是一個變量個數＊采樣個數的矩陣，這里3個信號，每個信號100列，則S為3*100
Sweight=rand（size（S1））;??%?取一隨機矩陣，作為信號混合的權矩陣
MixedS=Sweight*S;?????????%?得到三個信號的混合信號矩陣

%?將混合矩陣重新排列并輸出
subplot（334）plot（MixedS（1:））title（‘mixs1‘）axis（[0100-44]）;
subplot（335）plot（MixedS（2:））title（‘mixs2‘）axis（[0100-44]）;
subplot（336）plot（MixedS（3:））title（‘mixs3‘）axis（[0100-44]）;
MixedS_bak=MixedS;?????????????????????????%?將混合后的數據備份，以便在恢復時直接調用

%?標準化，預處理步驟之一
MixedS_mean=zeros（31）;
for?i=1:3
????MixedS_mean（i）=mean（MixedS（i:））;
end????????????????????????????????????????%?計算MixedS的均值
for?i=1:3
????for?j=1:size（MixedS2）
????????MixedS（ij）=MixedS（ij）-MixedS_mean（i）;
????end
end

%?白化，預處理步驟之一
MixedS_cov=cov（MixedS‘）;????????????????????%?cov為求協方差的函數
[ED]=eig（MixedS_cov）;??????????????????????%?對信號矩陣的協方差函數進行特征值分解
white=inv（sqrt（D））*（E）‘;????????????????????????%?Q為白化矩陣
MixedS_white=white*MixedS;??????????????????????%?MixedS_white為白化后的信號矩陣MixedS為混合信號矩陣
????
%?優化迭代，這里可以用多種方法
X=MixedS_white;????????????????????????????%?以下算法將對X進行操作X=MixedS_white為白化后的信號矩陣就是Z
[VariableNumSampleNum]=size（X）;
numofIC=VariableNum;???????????????????????%?在此應用中，獨立元個數等于變量個數信號個數，X=MixedS_white白化后信號矩陣的行數
B=zeros（numofICVariableNum）;?????????????%?初始化列向量w的寄存矩陣B=[b1??b2??...???bd]，初始化權值W矩陣，n*n，n為信號個數3個信號，3行3列
for?r=1:numofIC
????i=1;
????maxIterationsNum=100;???????????????%?設置最大迭代次數（即對于每個獨立分量而言迭代均不超過此次數）
????IterationsNum=0;
????b=rand（numofIC1）-.5;??????????????????%?隨機設置b初值3行1列
????b=b/norm（b）;???????????????????????????%?對b標準化?norm（b）:向量元素平方和開根號相當于單位化向量b
????while?i<=maxIterationsNum+1
????????if?i?==?maxIterationsNum???????????%?循環結束處理
????????????fprintf（‘\n第%d分量在%d次迭代內并不收斂。‘?rmaxIterationsNum）;
????????????break;
????????end
????????bb=b;??%?3*1???????????????????????
????????a2=1;
????????u=1;
????????t=X‘*b;?%?3*1
????????g=t.*exp（-a2*t.^2/2）;%相當于g=y.e（-y2/2）
????????dg=（1-a2*t.^2）.*exp（-a2*t.^2/2）;%對g的求導
????????b=（（1-u）*t‘*g*b+u*X*g）/SampleNum-mean（dg）*b;%對應公式3-11
???????????????????????????????????????????%?核心公式，參見理論部分公式2.52
????????b=b-B*B‘*b;????????????????????????%?對b正交化
????????b=b/norm（b）;?
????????if?abs（abs（b‘*bb）-1）<1e-19????????%?如果收斂，則
?????????????B（:r）=b;?????????????????????%?保存所得向量b
?????????????break;
?????????end
????????i=i+1;????????
????end
%????B（:r）=b;????????????????????????????????%?保存所得向量b
end

%%??信號重建，意思是將ICAedS與原始源信號S比較，會有差異，但是差異越小越好
ICAedS=B‘*white*MixedS_bak;??????%?計算ICA后的矩陣

%?將混合矩陣重新排列并輸出
subplot（33

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????3707??2013-09-24?13:03??fushangST1.m