wrls_go.m

来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 249 行

% *****************************************************************
%	wrls algorithm.(wrls_go.m)
%	Mind-Machine Interaction Research Center,
%	University of Florida.
%	
%	 Originally written  by:
%	John M. White, MSEE
%	10/28/91
%  Modified to fit the general speech analysis tool by Karthik
%  Date : July 18, 1997
%
%	NOTE: The effectiveness of the algorithm depends upon
%	the correct choice for the value of the variable "e1".
%	This is the value of the error, and is also denoted E1.
%	If e1 is too small, the VFF, lambda, will drop out excessively;
%	If e1 is too large, the VFF, lambda, will not drop out at all.
%	We presently choose a value for e1 manually, that is we try the
% 	algorithm several times with different values of e1 to find the
% 	correct order of magnitude.
%
% *****************************************************************
%
NEW_LEFT = 1;
NEW_RIGHT = length(SPEECH);

kk=[NEW_LEFT:NEW_RIGHT];
num_iter=NEW_RIGHT - NEW_LEFT+1;

theta_hat=zeros(vec_length, 1);
phi_hat=zeros(vec_length,1);

%Get values of variables from the window
lambda_min = eval(get(wrls_lambda_edit_h,'String'));
e1 = eval(get(wrls_e1_edit_h,'String'));



%%%%%%%%%%%%%%%%%%%%%%%%%
%window to display wait message
wait_window_1 = figure('Numbertitle','off',...
   'Color',BACK_COLOR,...
   'Name','Doing WRLS - Variable Forgetting Factor Analysis .. Please wait',...
   'Position',[117 400 560 70]);
pause(0.01);

wait_window_1_display = uicontrol('Style','Text',...
   'Position',[50 20 410 30],...
   'Backgroundcolor','white',...
   'Foregroundcolor','blue',...
   'String','Doing WRLS -VFF analysis. Please be patient');
pause(0.01);
%%%%%%%%%%%%%%%%%%%%%%%%%

for k=1:(wrls_num_poles-1),
	phi_hat(k,1)=-1.0*SPEECH(NEW_LEFT-1+wrls_num_poles-k);
end;

%Window to display calculation progress

wait_counter = 0;
wait_window = figure('Numbertitle','off',...
   'Color',BACK_COLOR,...
   'Name','Calculating',...
   'Position',[117 200 560 150]);
pause(0.01);

wait_window_slider = uicontrol('Style','Slider',...
   'Position',[50 100 410 30],...
   'Max',num_iter,...
   'Min',0,...
   'Value',wait_counter);

wait_window_display = uicontrol('Style','Text',...
   'Position',[50 50 410 30],...
   'Backgroundcolor','white',...
   'Foregroundcolor','blue',...
   'String','Calculating Frequency Response every tenth sample');



p=eye(vec_length);
p=p*1e6;

% ****************variables *************************
error=zeros(1,num_iter);           	% prediction error;
K=zeros(vec_length,1);    	% Gain;
lambda=zeros(1,num_iter);        	% forgetting factor;

% ****** initialize u, u_old ************************
u_hat=0;
u_hat_old=0;
u_hatt=zeros(1,num_iter);
surface=zeros(256,fix(num_iter/10));	%initialize 3D mesh values

% start main 

for k=(wrls_num_poles+1):(num_iter),
   u_hat_old=u_hat;
   theta_hat_old=theta_hat;
   p_old=p;
   %update phi_hat
   phi_hat((2:wrls_num_poles),1)=phi_hat((1:(wrls_num_poles-1)),1);
   phi_hat(1,1)=-1.0*SPEECH(NEW_LEFT-1+k-1);
   if wrls_num_zeros>=2,
      phi_hat((wrls_num_poles+2):(vec_length),1)=phi_hat((wrls_num_poles+1):(vec_length-1),1);
   end
   if wrls_num_zeros>=1,
      phi_hat((wrls_num_poles+1),1)=u_hat_old;
   end
   %update prediction error;
   error(k)=SPEECH(NEW_LEFT-1+k)-(phi_hat'*theta_hat_old);
   %update filter gain;
   K = p_old*phi_hat*inv(lambda(k-1)+phi_hat'*p_old*phi_hat);
   %update lambda;
   lambda(k)=1-(error(k)*error(k)*(1-phi_hat'*K)*(1-phi_hat'*K)/e1);
	if lambda(k)<lambda_min,
		lambda(k)=lambda_min;
	end
	%	update covariance;
	p=(1/lambda(k))*(p_old-(K*phi_hat'*p_old));
	%	input estimate;
	if lambda(k)<=(lambda_min+0.01),
		%pulse input;
		u_hat=SPEECH(NEW_LEFT-1+k)-phi_hat'*theta_hat_old;
	else
		%white input;
		u_hat=error(k)*(1-phi_hat'*K);
	end
	u_hatt(k)=u_hat;
	%	update parameter;
	theta_hat=theta_hat_old + K*error(k);
	
   %Display Progress
   wait_counter = k;
   set(wait_window_slider,'Value',wait_counter);
      
	%find frequency response;
   if rem(k,10)==0
      set(wait_window_display,'String','Calculating Frequency Response every tenth sample.');
      pause(0.01);
		[h,w]=freqz([1; theta_hat(wrls_num_poles+1:vec_length)],[1;theta_hat(1:wrls_num_poles)],256,1/Ts);
      surface(1:256,fix(k/10))=20*log10(abs(h(1:256,1)));

   end
end
set(wait_window_display,'String','Done with frequency response..Proceeding with main loop.');
set(wait_window_display,'String','Done with WRLS -VFF Analysis .. Now plotting the results.');
pause(0.01);

set(wait_window_1_display,'String','Done with WRLS -VFF Analysis .. Now plotting the results.');
pause(0.01);



close(wait_window);
clear wait_window;
clear wait_counter;
pause(0.01);


PV=[270 44 515 268];
s2 = 'WRLS - Waterfall';

% Open analysis window
while exist('wrls_out_win1_h')==1
   try1 = 'get(wrls_out_win1_h,''position'');';
   eval(try1,catch2);
   if check ==0
      clear wrls_out_win1_h;
      check = 1;
      break;
   end
   s1 = get(wrls_out_win1_h,'Name');
   if ~strcmp(s1,s2)
      clear wrls_out_win1_h;
      break;
   end
   figure(wrls_out_win1_h);
   break;
end;

if exist('wrls_out_win1_h')~=1;
	wrls_out_win1_h=figure('Position',PV,...
				'Resize','on',...
				'Numbertitle','off',...
				'Name',s2);
end
%Plotting Statements
vf=1:256;
vf=vf/(256)*1/(2*Ts);%frequency axis scaling
vh=1:fix(num_iter/10);%time axis scaling
vh=vh*10;
mesh(vf,vh,surface');%transpose surface to rearrange data
view(45,30);%rotate display to view it better
xlabel('Frequency');
ylabel('Time');
zlabel('Amplitude');
%mesh(surface,[-38 30]);
title('Log of Frequency Response');


PV=[270 290 515 275];
s2 = 'WRLS - Waveforms';
% Open analysis window
while exist('wrls_out_win2_h')==1
   try1 = 'get(wrls_out_win2_h,''position'');';
   eval(try1,catch2);
   if check ==0
      clear wrls_out_win2_h;
      check = 1;
      break;
   end
   s1 = get(wrls_out_win2_h,'Name');
   if ~strcmp(s1,s2)
      clear wrls_out_win2_h;
      break;
   end
   figure(wrls_out_win2_h);
   break;
end;

if exist('wrls_out_win2_h')~=1;
   wrls_out_win2_h=figure('Position',PV,...
      'Resize','on',...
      'Numbertitle','off',...
      'Name',s2);
end

subplot(411);
axis([NEW_LEFT NEW_RIGHT 0 1.2]);
plot(kk,lambda);
title('Lambda (Forgetting factor)');
subplot(412);
axis([NEW_LEFT NEW_RIGHT round(min(u_hatt)) round(max(u_hatt))]);
plot(kk,u_hatt);
title('Estimated input');
subplot(413);
axis([NEW_LEFT NEW_RIGHT round(min(error)) round(max(error))]);
plot(kk,error);
title('Error');
subplot(414);
axis([NEW_LEFT NEW_RIGHT round(min(SPEECH)) round(max(SPEECH))]);
plot(kk,SPEECH);
title('Speech');

close(wait_window_1);
clear wait_window_1;
pause(0.01);