独立成分分析 matlab 程序

FastICA算法的基本步骤：

1. 对观测辩唯数据进行中心化，使它的均值为0；

2. 对数据进游嫌行白化，。

3. 选择需要估计的分量的个数，设迭代次数

4. 选择一个初始权矢量（随机的）。

5. 令,非线性函数的选取见前文。

6. 。

7. 令。

8. 假如不收敛的话，返回第5步。

9．令，如果，返回第4步。

二．MATLAB源程序及说明：

%下程序为ICA的调用函数，输入为观察的信号，输出为解混后的信号

function Z=ICA(X)

%-----------去均值---------

[M,T] = size(X)%获取输入矩阵的行/列数，行数为观测数据的数目，列数为采样点数

average= mean(X')' %均值

for i=1:M

X(i,:)=X(i,:)-average(i)*ones(1,T)

end

%---------白化/球化------

Cx = cov(X',1) %计算协方差矩阵Cx

[eigvector,eigvalue] = eig(Cx)%计算Cx的特征值和特征向量

W=eigvalue^(-1/2)*eigvector' %白化矩阵

Z=W*X %正交矩阵

%----------迭代-------

Maxcount=10000 %最大迭代次数

Critical=0.00001 %判断是神灶手否收敛

m=M %需要估计的分量的个数

W=rand(m)

for n=1:m

WP=W(:,n) %初始权矢量（任意）

% Y=WP'*Z

% G=Y.^3%G为非线性函数，可取y^3等

% GG=3*Y.^2 %G的导数

count=0

LastWP=zeros(m,1)

W(:,n)=W(:,n)/norm(W(:,n))

while abs(WP-LastWP)&abs(WP+LastWP)>Critical

count=count+1 %迭代次数

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('未找到相应的信号)

return

end

end

W(:,n)=WP

end

Z=W'*Z

%以下为主程序，主要为原始信号的产生，观察信号和解混信号的作图

clear allclc

N=200n=1:N%N为采样点数

s1=2*sin(0.02*pi*n)%正弦信号

t=1:Ns2=2*square(100*t,50)%方波信号

a=linspace(1,-1,25)s3=2*[a,a,a,a,a,a,a,a]%锯齿信号

s4=rand(1,N)%随机噪声

S=[s1s2s3s4]%信号组成4*N

A=rand(4,4)

X=A*S%观察信号

%源信号波形图

figure(1)subplot(4,1,1)plot(s1)axis([0 N -5,5])title('源信号')

subplot(4,1,2)plot(s2)axis([0 N -5,5])

subplot(4,1,3)plot(s3)axis([0 N -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')

均衡器Matlab例程(1)解析

Equalizing Using a Training Sequence

例程代码

% Set up parameters and signals.

M = 4% Alphabet size for modulation

msg = randint(1500,1,M)% Random message，产生一个1500 x 1 在[0，M-1]区域内的随机整数序列

modmsg = pskmod(msg,M)% Modulate using QPSK. 进行QPSK调制的映射

trainlen = 500% Length of training sequence 定义训练序列的帧长

chan = [.986.845.237.123+.31i]% Channel coefficients 定义信道参量

filtmsg = filter(chan,1,modmsg)% Introduce channel distortion. 模拟信道变化

% Equalize the received signal.

eq1 = lineareq(8, lms(0.01))% Create an equalizer object.

eq1.SigConst = pskmod([0:M-1],M)% Set signal constellation. 设置星座图

[symbolest,yd] = equalize(eq1,filtmsg,modmsg(1:trainlen))% Equalize.

% Plot signals.

h = scatterplot(filtmsg,1,trainlen,'bx')hold on以蓝星画出未经过均衡的信号

scatterplot(symbolest,1,trainlen,'g.',h)在原图的基础上以绿色画出经过均衡后的信号

scatterplot(eq1.SigConst,1,0,'k*',h)在原图的基础上以黄色标出理想星座图

legend('Filtered signal','Equalized signal',...

'Ideal signal constellation')

hold off

% Compute error rates with and without equalization. 计算误码率

demodmsg_noeq = pskdemod(filtmsg,M)% Demodulate unequalized signal. 解调未均衡的码字

demodmsg = pskdemod(yd,M)% Demodulate detected signal from equalizer.解调已均衡的码字

[nnoeq,rnoeq] = symerr(demodmsg_noeq(trainlen+1:end),...对比未均衡信号和样本信号的误码和误码率

msg(trainlen+1:end))

[neq,req] = symerr(demodmsg(trainlen+1:end),...

msg(trainlen+1:end))对比均衡信号和样本信号的误码和误码率

disp('Symbol error rates with and without equalizer:')

disp([req rnoeq]

---