mgif_go.m

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

%
% Manuall glottal inverse filtering algorithm
% Author : Minkyu Lee
% Date 10-Oct-1994
% Modified by D. G. Childers
%

% get a differentiated glottal wave estimate

% do a fixed pre-emphasis
p_ydata_m=filter([1 -0.9], 1, ydata_m);
Fs_m=fft(p_ydata_m,1024);


LPF3=fir1(3,0.2);
ydata_lp=filter(LPF3,1,ydata_m);
p_ydata_lp=filter(LPF3,1,p_ydata_m);

% Calculate the transfer function from input ( Formants freq. and bandwidth )

if num_poles > 10
	ii=5;
else
	ii=ceil(num_poles/2);
end

for i=1:ii
	C(i)=-exp(-2*pi*fb(i)*Ts);
	B(i)=2*exp(-pi*fb(i)*Ts)*cos(2*pi*ff(i)*Ts);
	A(i)=1-B(i)-C(i);
end

nom=1;
for i=1:ii
	nom=nom*A(i);
end

den=conv([1 -B(1) -C(1)], [1 -B(2) -C(2)]);
for jj=3:ii
	den=conv(den,[1 -B(jj) -C(jj)]);
end

dg_m=filter(den, nom, ydata_lp);
dg_m=filter([1 -1],[1 -.99],dg_m);
if mgif_ps_flag == 0
	DGVV(NEW_LEFT:NEW_RIGHT)=dg_m(ext_left+1:length(dg_m)-ext_right);
elseif mgif_ps_flag == 1
	DGVV(gci_m(cur_gci)+NEW_LEFT:gci_m(cur_gci+1)+NEW_LEFT) = ...
      dg_m(ext_left+1:length(dg_m)-ext_right);
end

re_m=filter(den, nom, p_ydata_lp);
re_m=filter([1 -1],[1 -.99],re_m);
if mgif_ps_flag == 0
   RESIDUE(NEW_LEFT:NEW_RIGHT)=re_m(ext_left+1:length(re_m)-ext_right);
elseif mgif_ps_flag == 1
   RESIDUE(gci_m(cur_gci)+NEW_LEFT:gci_m(cur_gci+1)+NEW_LEFT) = ...
      re_m(ext_left+1:length(re_m)-ext_right);
end

% get a glottal wave estimate
g_m=zeros(size(dg_m));
isum=0;
for i=1:length(dg_m),
   isum=isum+dg_m(i);
   g_m(i)=isum;
end

g_m=filter([1 -1],[1 -.99],g_m);
if mgif_ps_flag == 0
   GVV(NEW_LEFT:NEW_RIGHT)=g_m(ext_left+1:length(g_m)-ext_right);
elseif mgif_ps_flag == 1
   GVV(gci_m(cur_gci)+NEW_LEFT:gci_m(cur_gci+1)+NEW_LEFT) = ...
      g_m(ext_left+1:length(g_m)-ext_right);
end

PV=[270 44 515 268];
s2 = 'MGIF - Output';

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

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

figure(mgif_out_win_h);
subplot(311),plot(RESIDUE);
title('Residue');
subplot(312),plot(DGVV);
title('Differentiated gvv (dgvv)');
subplot(313),plot(GVV);
title('Glottal volume velocity (gvv)');

figure(mgif_plot_win_h);
vh=1:1024/2;
vh=vh/(1024/2)*1/(2*Ts);
plot(vh, 20*log10(abs(Fs_m(1:512))));
hold on;
[H_m W_m]=freqz(nom, den,512);
plot(vh, 20*log10(abs(H_m))+gain,':');
p_dg_m=filter([1 -0.9], 1, dg_m);
F_m=fft(p_dg_m,1024);
plot(vh, 20*log10(abs(F_m(1:512)))-120,'-.');
hold off;

clc;