clear all
close all
clc
%展示
% [x,fs] = audioread('test01_16k.wav');
[x,fs] = audioread('recordta.wav');
t = [0:length(x)-1]/fs;
figure(1);
plot(t,x);
xlabel('t/s');
ylabel('amp');
title('时域波形');
%清音浊音片段
t1 = 0.27;%清音开始
t2 = 0.34;%清音结束
t3 = 0.40;%浊音开始
t4 = 0.60;%浊音结束
%%%%%%%%%%同态滤波部分
unvoice_x = x([ fix(t1*fs) :fix( t2*fs)]);
voice_x = x([ fix(t3*fs) :fix( t4*fs)]);
x = voice_x;
index1 =1;
nfft = 1024;
index2=index1+nfft-1;
if( index2 > length(x))
error('信号太短')
end
x=x(index1:index2);
figure(2)
plot(x);
%倒谱
z=rceps(x); %MATLAB提供的倒频谱函数
figure(3)
plot(z(1:512));
%选择倒谱提升窗的边界点
mcep=40;
%声道信息
zz1 = zeros(size(z));
%激励信息
zz2 = zeros(size(z));
zz1(1:mcep) = z (1:mcep);
zz1(end - mcep +1 :end) = z (end - mcep +1 :end);
zz2(mcep+1:511) = z (mcep+1:511);
zz2(512:end-mcep) = z (512:end-mcep);
%声道信息
% z1=real(ifft(exp(fft(zz1))));
z1= fft(real(ifft(exp(fft(zz1)))));
%激励信息
z2=real(ifft(exp(fft(zz2))));
figure
subplot(211)
plot(abs(z2(1:512)));
title('激励信号');
xlabel('时域');
ylabel('幅度');
% ylim([0 5]);
subplot(212)
plot(20*log10(abs(z1(1:512))));
title('声门信号');
xlabel('频域');
ylabel('分贝');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%LPC参数解卷积
% 得到元音（浊音）
% x = voice_x;
x = voice_x;
% 最小频率点50hz
ms20=fs/50;
% 计算自相关
r=xcorr(x,ms20,'coeff');
%最大频率点
ms2=fs/500 ;
%最小频率点
ms20=fs/50 ;
%确定delay
r=r(ms20+1:2*ms20+1);
[rmax,tx]=max(r(ms2:ms20));
f0 = fs/(ms2+tx-1);
flag = 0;
%相关性太小，认为不是周期信号
if( rmax < 0.6)
A = (rmax);
fprintf('清音信号\n');
flag = 0;
else
A = (rmax);
fprintf('浊音信号\n');
fprintf('最大值rmax=%g 基频Fx=%gHz\n',rmax,f0);
flag =1;
end
%
%得到线性预测LP系数
ncoeff=2+fs/1000;
%lpc估计
a=lpc(x,ncoeff);
% [h t]= impz(1,a);
% figure
% plot(t,h);
%打印频谱
[h,f]=freqz(1,a,512,fs);
figure
subplot(211)
t = [0:length(x)-1]/fs;
% sim = A*sin(2*pi*f0*t);
f0 = fs/fix(fs/f0);
%creating pulsetrain;
for i =1:1024
if i./((ms2+tx-1)) == floor(i./((ms2+tx-1)))
sim(i) = 1;
else
sim (i) = 0;
end
end
plot(sim)
ylim([-2,2])
xlabel('时间');
ylabel('激励信号amp');
if (flag==0)
title('清音,激励信号不做参考');
else
title('浊音,激励信号');
end
% suptitle('参数提取方法');
subplot(2,1,2);
plot(f,20*log10(abs(h)+eps));
% plot(f,abs(h));
legend('LPC函数');
xlabel('频域');
ylabel('增益(dB)');
if (flag==0)
title('清音，参数LPC声道不做参考');
else
title('浊音，LPC声道信息');
end
% ylim([-1 1])
% NFFT = length(sim);
% Y = abs(fft(sim,NFFT))/NFFT;
% f = fs/2*linspace(0,1,NFFT/2+1);
