matlab语音信号盲分离.doc

武汉理工大学《信息处理课群综合训练与设计》说明书

目录
摘要 1
abstract 2
1 ICA 的基本原理及特点 3
1.1 盲分离数学模型 3
1.2 ICA 算法描述 4
1.3 FICA算法 4
2 FICA设计思想 5
3 实验仿真结果记录 6
3.1 仿真时域波形及频谱 6
3.1.1 原始信号 6
3.1.2 混合信号 9
3.1.3 分离信号 12
3.2 仿真所用的源程序 15
4 实验结果分析 24
5 小结与体会 25
6 参考文献 26

摘要
语音信号的盲源就是在源信号和源信号如何混合都未知的情况下,从观测到的混合信号中恢复出未知源信号。语音信号盲分离技术被成功地用在了通信、医学、图像和语音信号处理等领域。我们所要研究的混合语音信号盲分离问题就是用麦克风阵列或多个麦克风阵列来模仿人的耳朵,采集得到相互干扰的混叠语音信号,然后通过分离算法将混叠的语音信号相互分离开来,提取我们所感兴趣的信号。举个例子就是在多人同时说话的嘈杂环境下,我们能够辨识感兴趣人的说话声的能力。然后把它分辨出来。

abstract
The speech signal blind source is in the source signal and the source signal how to mix all unknown, from observations of mixed signal in recovering the unknown source signal. Blind separation technology of speech signal was successfully used in the communication, medical, image and voice signal processing, etc. We have to study mix of speech signal is blind source separation with a microphone array or more microphone array to imitate human ears, acquisition get each other interference aliasing speech signal, and then through the separation algorithm will aliasing voice signal mutual separated, extraction we are interested in signal. For example is more than in noisy environment of speak at the same time, we can identify interested in the ability of the human voice. And then take it apart.

1 ICA 的基本原理及特点
1.1 盲分离数学模型
盲信号分离是指在没有任何先验知识的条件下，仅根据源信号之间的统计独立特性和由传感器输出的观测信号，把源信号分离出来。独立分量分析正是盲信号分离的有效方法之一。
我们将该问题表述成图一所示的线性模型
图一 盲信号分离的数学模
该线性模型可以用下式表示：
x = A⋅s..............................(公式一)
u = W⋅x..............................(公式二)
源信号s=(s1,s2......sn)和混合矩阵A都是未知的，W为待求分离矩阵，经过分离矩阵W 的作用，最终得到分离后的输出矢量u=(u1,u2......un);观测信号x=(x1,x2.......xn)已知。盲信号分离的最终目的是为了寻找分离矩阵W，使输出信号u尽可能地逼近源信号s。
以两个信号源为例用恒等式表示：
x1(t)= h11s(t) + h12j(t).......................(公式三)
x2(t)= h21s(t) + h22j(t)......................（公式四）
y(t)= w1x1(t) + w2x2(t).......................（公式五）
其中x1(t) ，x2 (t) 表示两个未知混合信号，s(t) ，j(t)表示收集到的两个信号源，h11 ， h12 ，h21 ，h22表示混合矩阵，不同个数的信号源输入对应不同阶数的矩阵。而y(t) 则表示最后解混后的声音，w1 ，w2 表示分离矩阵，我们得主要任
务选取合适的分离矩阵，使得最后的输出y(t)最大限度的接近s(t) 或者j(t)，亦既是我们所需要挺清楚的没有其它干扰的信号声源。
1.2 ICA 算法描述
如果我们依据上面提出的模型，对分离过程和分离结果没有任何限定和约束，那么显然解并不是唯一的；而我们却希望输出尽量接近于所期望的结果。独立分量分析(ICA)为这类问题提供了新的解决思路。
ICA的基本思想是假设源信号彼此相互独立，寻找分离矩阵的目的是使输出结果的各分量尽可能地相互独立，我们可以把它看作是一个最优化问题，目标函数用来度量分离结果的独立性。
从处理技术上看，依据独立性作分解势必涉及概率密度函数和高阶统计量；而且处理过程中常常要引入非线性环节。从这一意义上看，ICA算法思想优于常用的只建立在二阶统计量上的线性处理技术。
进行ICA分解的基本原则可以简单地概括为两条：
（1）非线性去相关：求分离矩阵使其任意两输出不但本身不相关，而且经非线性变换后的分量也不相关。这里的目标函数的选择对分离效果至关重要。
（2）使输出尽可能非高斯化：在输出某分量的方差一定的条件下，将输入各分量作线 性组合。优化选择各权重，使输出分量尽可能非高斯化。输出的非高斯性的每一个局部最大值即给出一个独立分量。
1.3 FICA算法
与普通 ICA 相比较而言，FICA 具有更快的收敛速度，而且还省去了选择步长常数的麻烦，迭代次数也相应减少。尤其是在处理分离复杂多变的语音信号时，优势非常明显。
FICA对观测信号x进行预处理，包括两部分:去均值化和预白化，即通过对观测信号进 行线性变换，得到一个一个新的白噪声化矢量z，使其各个分量不关联。这样可以大大简化ICA算法。所以此次课程设计我采用FICA算对三路语音信号进行盲处理。

2 FICA设计思想
我们知道语音信号是非平稳的信号，但是它在团时间内是具有稳定性的，我根据根据盲信号分离原理，用matlab采集三路语音信号，分别为music.wav,dragen.wav以及man.wav。选择合适的混合矩阵生成三路混合信号。利用盲信号分离算法中的独立成分分析ICA进行语音盲信号的分离，将源信号分别从混合信号中分离出来。
按照盲分离模型FICA具体步骤如下：
（1）将三段声音文件保存为“.wav”格式，分别以列向量形式读入，记录时域波形以及频谱图。
（2）对源信号进行预处理，方差为1，均值为0，归一化处理后得到归一化频谱。
（3）根据X = A⋅S，混合矩阵A选取二阶随机矩阵，A=rand(2),经A作用后，生成混合语音序列X。
（4）对X进行方差1，均值为0的操作，得到归一化频谱；然后进行预白化处理，得到白化后的信号。

3 实验仿真结果记录
3.1 仿真时域波形及频谱
3.1.1 原始信号

图二 man.wav原始时域波形图

图三 man.wav原始频谱图

图四 dragen.wav原始时域波形图

图五 dragen.mav原始频谱图

图六 music.wav原始时域波形图

图七 music.wav原始频谱图
3.1.2 混合信号

图八 混合语音信号man时域波形图

图九 混合语音信号man频谱图

图十 混合语音信号dragen时域波形图

图十一 混合语音信号dragen频谱图

图十二 混合语音信号music时域波形图

图十三 混合语音信号music频谱图
3.1.3 分离信号

图十四 分离语音信号man时域波形图

图十五 分离语音信号man频谱图

图十六 分离语音信号dragen时域波形图

图十七 分离语音信号dragen频谱图

图十八 分离语音信号music时域波形图

图十九 分离语音信号music频谱图
3.2 仿真所用的源程序
% data length
N = 400000;
% read data from wave file
s1= wavread('C:\Users\Administrator\Desktop\临时\37724126FastICA\Fast ICA\music.wav');
s1=s1(1:N)';
s2=wavread('C:\Users\Administrator\Desktop\临时\37724126FastICA\Fast ICA\man.wav');
s2=s2(1:N)';
s3=wavread('C:\Users\Administrator\Desktop\临时\37724126FastICA\Fast ICA\dragen.wav');
s3=s3(1:N)';
s=[s1; s2; s3];

% mix the data
A=[1 4 5
3.4 2.4 2
6 2 3.5];
x = A * s;

%recover s1, s2 & s3 from x
len = size(x, 2); % number of data samples
IC_num = size(x, 1); % number of ICs

% Center the data with 0 mean
for cnt = 1: IC_num
x(cnt, :) = x(cnt, :) - mean(x(cnt, :)')';
end

% Whiten the data
Rx = x*x'/len;
[Ux Sx] = eig(Rx); % assert that Ux*Sx*Ux'-Rx = 0

Vx = inv(Sx).^0.5*Ux'; % Vx is the whiten matrix
z = Vx * x; % z is whiten data
% assert that z*z'/len = I
% Recover ICs
w = eye(IC_num);
w_old = w;

eps = 1.0e-5;
flag = 1;
while (flag == 1)
for (cnt = 1: IC_num)
y = w(:, cnt)'*z;
Gy = -exp(-y.^2/2);
gy = y.*(-Gy);
gy1 = (1-y.^2).*(-Gy);
E1 = z*gy'/len;
E2 = mean(gy1')';
w(:, cnt) = E1-E2*w(:, cnt);
end
wt = w';
[Uw Sw] = svd(wt*wt');
wt = inv(Uw * sqrt(Sw)) * wt;
w = wt';
change=((sum(sum((abs(w'*w_old) - eye(IC_num)).^2)))/IC_num);
if ( change < eps)
flag = 0;
end
w_old = w;
End
Display(w)

% Recover s
s_n = w'*z;

% Play all the music
%play Original sound 1
[x1,fs,bits]=wavread('C:\Users\Administrator\Desktop\临时\37724126FastICA\Fast ICA\man.wav'); %读取语音信号的数据，赋给变量x1
disp('press a key to hear Original sound 1');
pause;
sound(x1,fs); %播放语音信号

%Show time domain waveform figure and spectrum image 1
y1=fft(x1,1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(x1) %做原始语音信号的时域图形
title('原始语音信号man波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图

plot(f,abs(y1(1:512)));
title('原始语音信号man频谱')
xlabel('Hz');
ylabel('幅值');

%play Original sound 2
[x2,fs,bits]=wavread('C:\Users\Administrator\Desktop\临时\37724126FastICA\Fast ICA\dragen.wav'); %读取语音信号的数据，赋给变量x1
disp('press a key to hear Original sound 2');
pause;
sound(x2,fs); %播放语音信号
%Show time domain waveform figure and spectrum image 2
y2=fft(x2,1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(x2) %做原始语音信号的时域图形
title('原始语音信号dragen波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图
plot(f,abs(y2(1:512)));
title('原始语音信号dragen频谱')
xlabel('Hz');
ylabel('幅值');

%play Original sound 3
[x3,fs,bits]=wavread('C:\Users\Administrator\Desktop\临时\37724126FastICA\Fast ICA\music.wav'); %读取语音信号的数据，赋给变量x1
disp('press a key to hear Original sound 3');
pause;
sound(x3,fs); %播放语音信号
%Show time domain waveform figure and spectrum image 3
y3=fft(x3,1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(x3) %做原始语音信号的时域图形
title('原始语音信号music波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图
plot(f,abs(y3(1:512)));
title('原始语音信号music频谱')
xlabel('Hz');
ylabel('幅值');

% Play mix sound 1
x(1,:)=x(1,:)/(max(s_n(1,:))-min(x(1,:)))*2;
disp('press a key to hear mix sound 1');
pause;
wavplay(x(1, :), 44100);
%Show time domain waveform figure and spectrum image 4

y4=fft(x(1,:),1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(x(1,:)) %做原始语音信号的时域图形
title('混合语音信号man波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图
plot(f,abs(y4(1:512)));
title('混合语音信号man频谱')
xlabel('Hz');
ylabel('幅值');

% Play mix sound 2
x(2,:)=x(2,:)/(max(x(2,:))-min(x(2,:)))*2;
disp('press a key to hear mix sound 2');
pause;
wavplay(x(2, :),44100);

%Show time domain waveform figure and spectrum image 5
y5=fft(x(2,:),1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(x(2,:)) %做原始语音信号的时域图形
title('混合语音信号dragen波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图
plot(f,abs(y5(1:512)));
title('混合语音信号dragen频谱')
xlabel('Hz');
ylabel('幅值');

% Play mix sound 3
x(3,:)=x(3,:)/(max(x(3,:))-min(x(3,:)))*2;
disp('press a key to hear mix sound 3');
pause;
wavplay(x(3, :),44100);

%Show time domain waveform figure and spectrum image 6
y6=fft(x1,1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(x(3,:)) %做原始语音信号的时域图形
title('混合语音信号music波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图
plot(f,abs(y6(1:512)));
title('混合语音信号music频谱')
xlabel('Hz');
ylabel('幅值');
% Play sound 1
s_n(1,:)=s_n(1,:)-mean(s_n(1,:));
s_n(1,:)=s_n(1,:)/(max(s_n(1,:))-min(s_n(1,:)))*2;
disp('press a key to hear sound 1');
pause;
wavplay(s_n(1, :), 44100);

%Show time domain waveform figure and spectrum image 7
y7=fft(s_n(1,:),1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(s_n(1,:)) %做原始语音信号的时域图形
title('分离语音信号man波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图
plot(f,abs(y7(1:512)));
title('分离语音信号man频谱')
xlabel('Hz');
ylabel('幅值');

% Play sound 2
s_n(2,:)=s_n(2,:)-mean(s_n(2,:));
s_n(2,:)=s_n(2,:)/(max(s_n(2,:))-min(s_n(2,:)))*2;
disp('press a key to hear sound 2');
pause;
wavplay(s_n(2, :), 44100);

%Show time domain waveform figure and spectrum image 8
y8=fft(s_n(2,:),1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(s_n(2,:)) %做原始语音信号的时域图形
title('分离语音信号gragen波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图
plot(f,abs(y8(1:512)));
title('分离语音信号dragen频谱')
xlabel('Hz');
ylabel('幅值');

% Play sound 3
s_n(3,:)=s_n(3,:)-mean(s_n(3,:)); s_n(3,:)=s_n(3,:)/(max(s_n(3,:))-min(s_n(3,:)))*2; disp('press a key to hear sound 3');
pause;
wavplay(s_n(3, :), 44100);

%Show time domain waveform figure and spectrum image 9
y9=fft(s_n(3,:),1024); %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(s_n(3,:)) %做原始语音信号的时域图形
title('分离语音信号music波形');
xlabel('时间 n');
ylabel('幅值 n');
figure(2)
%做原始语音信号的FFT频谱图
plot(f,abs(y9(1:512)));
title('分离语音信号music频谱')
xlabel('Hz');
ylabel('幅值');

4 实验结果分析
由实验最终结果听到的声音以及波形图、频谱图的结果可以看出，最终分离出来的声音波形以及频谱图虽然有一定的区别，但是大体的分布还是一样的，说明利用盲信号分离算法中的独立成分分析ICA进行语音盲信号的分离，可以将源信号分别从混合信号中分离出来，即可以减弱不想要听到的声音，刷选出我们想要的那个声源。
这次设计我所选的混合矩阵为：
A=
[ 1 4 5
3.4 2.4 2
6 2 3.5 ]

求的的分离矩阵为：
W =
[ 0.9806 -0.1945 -0.0250
-0.0948 -0.3587 -0.9286
0.1717 0.9130 -0.3702]

5 小结与体会
这次课程设计我通过学习FastICA算法的理论知识，我初步的掌握的FastICA算法的原理和基本步骤。并且通过仿真，实现了语音信号和图象信号的分离。发现用FastICA算法实现语音和图象信号分离，结果还是比较理想的。
同时也对知识的学习有所感触。在大学期间，接触Matlab软件也就是在做信号与线性系统实验和数字信号处理实验中用到，当初并没有什么深入的学习。通过这次课程设计设计，我了解了Matlab语言，发现其要比C语言相对简单多。但是Matlab语言中的函数还是比较多的，在编程的时候，一定要借助Matlab帮助。另外当一个函数看不懂时，要把它代到程序中去，单独实现结果，通过结果来分析函数，象像很多函数是单凭看理论是看不懂的。另外，当看懂一个函数时，一定要在程序中注释，否则过一段时间就忘记了其功能。

6 参考文献
[1] 杨行峻，郑君里；《人工神经网络与盲信号处理》[M]. 北京: 清华大学出版社，2003.
[2] 王涛，周荷琴，冯焕清.《基于独立分量分析的混合声音信号分离》.中国科学技术大学学报， 2001 Vol.31 No.01   
[3] 赵浩,周卫东,钟凌惠.《独立分量分析在生物医学信号处理中的应用》.生物医学工程研究， 2003 Vol.22 No.4
[4] 童峰,陆佶人.《基于盲源分离的水声信道盲均衡处理方法》.电子与信息学报,2003 Vol.25 No.2
[5] 张贤达.《盲信号处理几个关键问题的研究.深圳大学学报（理工版）》， 2004 Vol.21 No.3
[6] 刘琚,何振亚.《盲源分离和盲反卷积》.电子学报， 2002 Vol.30 No.4
[7] 杨福生.《独立分量与分析与应用》. 清华大学出版社2005。

26