数字信号处理课程项目

1、数字信号处理B课程项目实施报告题 目： 数字音效处理器 . 目录一、设计要求：3二、盲信号处理3（1）Fast ICA算法简介3（2）遇到问题5三、去噪5（1）遇到问题5（2）维纳滤波器6四、性能与结果分析71维纳滤波7（1）去噪时域图7（2） 去噪频域图8（3） 频域瀑布图82盲信号分离9（1）分离男声9（2） 分离女声93用户界面10五、源程序代码11（1）FastICA盲信号分离代码：11（2）维纳滤波去噪：13（3） 瀑布图代码：15（4） 主界面代码：16七、小组分工17一、设计要求：设计一个基于MATLAB的数字音频分析与处理系统，能够实现对数字音频的测试分析与处理。要求：（1）

2、输入语音信号源为实际环境采集语音；（2） 至少实现1种音效测试分析功能（频率响应，瀑布频谱图，相位响应曲线，抗阻曲线，谐波失真曲线，互调失真曲线，音压曲线等）；并至少实现1中语音处理功能（声纹识别，关键词识别，语音去噪，声源分离等）；并分别对其进行实际运用。（3） 用人机交互界面操控、扬声器/耳机输出音效。二、盲信号处理（1）Fast ICA算法简介快速ICA算法，由于采用了定点迭代的思想，所以又称固定点(Fixed-Point)算法或定点算法，是一种快速寻优神经网络算法。与普通的神经网络算法不同的是这种算法采用了牛顿迭代的思想对采样点采用批处理的方式，在每一步迭代中有大量的样本数据参与运算。

3、基于目标函数的选择不同，从而产生了FastICA算法的基于负熵最大、峭度、互信息最小、最大似然估计等的多种形式。FastICA算法的优点如下：1.收敛速度快，FastICA收敛速度为2次以上，普通的ICA算法收敛仅为1次。2.能利用任何的非线性函数直接找出任何非高斯分布的独立分量，而对于其他的算法来说，概率密度函数的估计不得不首先进行。3.独立分量可被逐个估计出来，类似于做投影追踪，这在仅需要估计几个（不是全部）独立分量的情况下，能减小计算量。4.Fast-ICA算法有许多神经算法的优点，它是并行、分布式的且计算简单，内存要求很少，它的性能能够通过选择适当的非线性函数来最佳化。从Fast-IC

4、A算法流程图可以看出，通常先对数据进行初步的白化或球化处理，白化处理可去除各观测信号之间的相关性，从而简化了后续独立分量的提取过程。这是由于一般情况下，所获得的数据都具有相关性。通过去均值和白化等预处理在对其进行真正的算法分离。通常情况下，数据进行白化处理与不对数据进行白化处理相比，算法的收敛性比较好。概括来说，Fast-ICA算法主要包含两部分:预处理部分和核心算法部分。预处理部分主要包括去均值化和白化处理；核心算法部分是基于拟牛顿算法，具体的实现见附录程序。 Fast-ICA算法流程图:观测信号x去均值对去均值后的观测信号x进行白化基于负熵最大化理论确定目标函数应用牛顿迭代算法求

5、目标函数的最大值，即分离矩阵w收敛？由获得信号源是否在所有等方差的随机变量中，高斯变量的熵最大，因而我们可以利用熵来度量非高斯性，常用熵的修正形式，即负熵。在分离过程中，可通过对分离结果的非高斯性度量来表示分离结果间的相互独立性，当非高斯性度量达到最大时，则表明已完成对各独立分量的分离。负熵的定义： ，H为随机变量的微分熵当Y具有高斯分布时Ng（Y）=0，Y的非高斯性越强，其微分熵越小，所以Ng（Y）可以作为随机变量Y非高斯性的测度。快速ICA学习规则是找一个方向以便具有最大的非高斯性。FastICA算法的迭代公式：（2）遇到问题对于声源分离的理解一开始我们将其理解为经典的“鸡尾酒会效应”，经

6、过大量的数据查询和整理，我们发现目前情况下，我们只能用机器实现分离简单叠加的盲信号，或者说是以特定形式混合在一起的信号，不能分离一些复杂的盲信号。在近几年的有效方法中，FastICA算法可以有效的分离简单叠加的盲信号，并且可以用matlab编程的方式得以实现。于是，我们使用matlab中的FastICA算法进行了简单叠加的盲信号分离。三、去噪（1）遇到问题在本次课程项目中，我们一开始选用了最基础的低通滤波器，并花了一些时间做出了巴特沃斯低通滤波器进行测试。但是除了语音变得更沉闷一些了对噪音并没有什么实质性的消除效果。因为我们组要消除的噪音是雪花噪音。我们尝试了维纳滤波器，维纳滤波器可以将大部分

7、噪音都滤除，只留下需要的人声信号以及几乎听不见的混在人声特定地方的噪音。（2）维纳滤波器维纳（Wiener）是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看作是一种估计问题或一种线性估计问题。 一个线性系统，如果它的单位样本响应为，当输入一个随机信号，且 （1）其中表示信号，)表示噪声，则输出为 （2）我们希望通过线性系统后得到的尽量接近于，因此称为的估计值，用表示，即 （3） 则维纳滤波器的输入输出关系可用下面图1表示。图1实际上，式（2）所示的卷积形式可以理解为从当前和过去的观察值，来估计信号的当前值。因此，用进行过滤问题实际上是一种统计估计问题。一般地，从当前

8、的和过去的观察值，估计当前的信号值成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值称为外推或预测；从过去的观察值，估计过去的信号值称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准则的。如果我们分别以与表示信号的真实值与估计值，而用表示他们之间的误差，即 （4）显然可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小： （5）四、性能与结果分析1维纳滤波（1）去噪时域图从无语音段提取出的噪声波形维纳滤波后波形 信噪比=14.8306dB带噪语音

9、 信噪比=4.8483dB（2） 去噪频域图带噪语音频谱噪声频谱处理后的语音频谱（3） 频域瀑布图2盲信号分离（1）分离男声时域图频域图（2） 分离女声时域图频域图3用户界面五、源程序代码（1）FastICA盲信号分离代码：x,fs=audioread('D:shuzixinhaogongchengXMhunhe1.wav');y,fs=audioread('D:shuzixinhaogongchengXMhunhe2.wav');%sound(x,fs);  % /播放原语音信号N1=length(x);  %

10、 /求出语音信号的长度N2=length(y);t1=(0:N1-1)/fs;    % /语音信号的长度除以采样频率t2=(0:N2-1)/fs;w1=fs*linspace(0,1,N1); w2=fs*linspace(0,1,N2); x1=abs(fft(x,N1);  y1=abs(fft(y,N2);  x=x(:,1);y=y(:,1);x=x'y=y'MixedS = x;y;MixedS_bak=MixedS; 

11、;  % 将混合后的数据备份，以便在恢复时直接调用%  中心化  %MixedS_mean=zeros(2,1);for i=1:2    MixedS_mean(i)=mean(MixedS(i,:);endfor i=1:2    for j=1:size(MixedS,2)        MixedS(i,j)=MixedS(i,j

12、)-MixedS_mean(i);    endend%  白化  %MixedS_cov=cov(MixedS');           % cov为求协方差的函数E,D=eig(MixedS_cov);            % 对信号矩阵的协方差函数进行特

13、征值分解Q=inv(sqrt(D)*(E)'           % Q为白化矩阵MixedS_white=Q*MixedS;     % MixedS_white为白化后的信号矩阵IsI=cov(MixedS_white');   % IsI应为单位阵  %FASTICA算法  %X=MixedS_white; 

14、;        % 以下算法将对X进行操作VariableNum,SampleNum=size(X);numofIC=VariableNum;   % 在此应用中，独立元个数等于变量个数B=zeros(numofIC,VariableNum);  % 初始化列向量w的寄存矩阵,B=b1  b2  .   bdfor r=1:numofIC   &

15、#160;i=1;maxIterationsNum=100;                % 设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数）    b=rand(numofIC,1)-.5;                &

16、#160;  % 随机设置b初值    b=b/norm(b);                           % 对b标准化 norm(b):向量元素平方和开根号    while i<

17、=maxIterationsNum+1        bOld=b;                                  t=X'*b;  

18、;      g=t.*exp(-t.2/2);        dg=(1-t.2).*exp(-t.2/2);        b=X*g/SampleNum-mean(dg)*b;        b=b-B*B'*b;      

19、           % 对b正交化        b=b/norm(b);         if abs(abs(b'*bOld)-1)<1e-9    % 如果收敛，则     &#

20、160;       B(:,r)=b;          % 保存所得向量b             break;         end     

21、   i=i+1;            endendif i = maxIterationsNum+1          % 循环结束处理      fprintf('n第%d分量在%d次迭代内并不收敛。', r,maxI

22、terationsNum);      return;end%  ICA计算的数据复原并构图  %ICAedS=B'*Q*MixedS_bak;          % 计算ICA后的矩阵% 将混合矩阵重新排列并输出sound(ICAedS(1,:),fs);Ws=sum(ICAedS(1,:).2);%计算信噪比Wn=sum(ICAedS(2,:).2);A=B'

23、;*Q;a11=A(1,1);a12=A(1,2);a21=A(2,1);a22=A(2,2);SNRy1=10*(a112*Ws)/(a122*Wn)（2）维纳滤波去噪：if (nargin<3 | isstruct(IS)   % 如果输入参数小于3个或IS是结构数据   IS=.25; endW=fix(.025*fs);             

24、;  % 帧长为25msSP=.4;                                      % 帧移比例取40%(10ms) wnd=hamming(

25、W);                % 设置窗函数% 如果输入参数大于或等于3个并IS是结构数据(为了兼容其他程序)if (nargin>=3 && isstruct(IS)    SP=IS.shiftsize/W;    nfft=IS.nfft;   

26、; wnd=IS.window;    if isfield(IS,'IS')        IS=IS.IS     else        IS=.25;    endendpre_emph=0;signal=filter(1 -pre_emph,1,sign

27、al);      % 预加重NIS=fix(IS*fs-W)/(SP*W) +1);           % 计算无话段帧数y=segment(signal,W,SP,wnd);            % 分帧 Y=fft(y);     

28、;                              % FFTYPhase=angle(Y(1:fix(end/2)+1,:);         % 带噪语音的相位角 Y=abs(

29、Y(1:fix(end/2)+1,:);            % 取正频率谱值numberOfFrames=size(Y,2);             % 计算总帧数FreqResol=size(Y,1);          

30、     % 计算频谱中的谱线数N=mean(Y(:,1:NIS)')'                % 计算无话段噪声平均谱值 LambdaD=mean(Y(:,1:NIS)').2)'          % 初始噪声功率谱方差alpha=.99;

31、60;                    % 设置平滑系数fn=numberOfFrames;miniL=5;                        

32、;            % 设置miniLvoiceseg,vosl,SF,Ef=pitch_vad1(y,fn,T1,miniL); %端点检测 NoiseCounter=0;                 % 初始化NoiseCounterNoiseLength

33、=9;                  % 设置噪声平滑区间长度G=ones(size(N);                % 初始化谱估计器Gamma=G;X=zeros(size(Y);    

34、;            % 初始化Xh=waitbar(0,'Wait.');            % 设置运行进度条图 for i=1:numberOfFrames    SpeechFlag=SF(i);    if

35、 i<=NIS                               % 若i<=NIS在前导无声(噪声)段        SpeechFlag=0;  

36、;      NoiseCounter=100;    %else                   % i>NIS判断是否为有话帧        %NoiseFlag, SpeechFlag,

37、0;NoiseCounter, Dist=vad(Y(:,i),N,NoiseCounter);     end    if SpeechFlag=0              % 在无话段中平滑更新噪声谱值        N=(NoiseLength*N+Y

38、(:,i)/(NoiseLength+1);         LambdaD=(NoiseLength*LambdaD+(Y(:,i).2)./(1+NoiseLength);%更新和平滑噪声方差    end     gammaNew=(Y(:,i).2)./LambdaD;          % 计算后验

39、信噪比    xi=alpha*(G.2).*Gamma+(1-alpha).*max(gammaNew-1,0); % 计算先验信噪比    Gamma=gammaNew;    G=(xi./(xi+1);         % 计算维纳滤波器的谱估计器    X(:,i)=G.*Y(:,i);  

40、0;       % 维纳滤波后的幅值%显示运行进度条图    waitbar(i/numberOfFrames,h,num2str(fix(100*i/numberOfFrames);endclose(h);                       &#

41、160;           % 关闭运行进度条图output=OverlapAdd2(X,YPhase,W,SP*W);        % 语音合成output=filter(1,1 -pre_emph,output);      % 消除预加重影响output=output/max(abs(output)

42、; （3） 瀑布图代码：b,Fs=audioread('D:shuzixinhaogongchengXMxuehua.wav');b=b(:,1);t=size(b)/44100;f=0:1:22049;R=1;G=0;B=0.5;figurefor i=1:0.2:t/2; z=b(i*44100-44099):(i*44100); z=abs(fft(z,44100)/22050; z=z' f,i=meshgrid(f,i); plot3(f,i,z(1:22050),'color',R G B); hold on; G=G+0.035;endR=0.5;G=1;B=0;for i=t/2:0