【杨勇】小论文MATLAB实验程序主要部分代码--语音增强--stsa_weuclid1
2016-06-16 00:00
721 查看
摘要: stsa_weuclid1
function stsa_weuclid1(filename,outfile,p)
if nargin<3
fprintf('Usage: stsa_weuclid inFile outFile.wav p \n\n');
return;
end
fp=fopen(filename,'r');
if fp<=0, error('Unable to open input file'); end;
% ------------ get the file extension, and then read the whole file -----------
%
ind1=find(filename == '.');
if length(ind1)>1, ind=ind1(length(ind1)); else, ind=ind1; end;
ext = lower(filename(ind+1:length(filename)));
[HDRSIZE, Srate, bpsa, ftype] = gethdr(fp,ext);
x=fread(fp,inf,'short');
fclose(fp);
%x=x/32768;
% =============== Initialize variables ===============
%
len=floor(20Srate/1000); % Frame size in samples
if rem(len,2)==1, len=len+1; end;
PERC=50; % window overlap in percent of frame size
len1=floor(lenPERC/100);
len2=len-len1;
win=hanning(len); %tukey(len,PERC); % define window
U=norm(win);
% Noise magnitude calculations - assuming that the first 6 frames is noise/silence
%
nFFT=len;
nFFT2=len/2;
noise_mean=zeros(nFFT,1);
j=1;
for k=1:5
noise_mean=noise_mean+abs(fft(win.*x(j:j+len-1),nFFT)/U);
j=j+len;
end
noise_mu=noise_mean/5;
noise_mu2=noise_mu.^2;
%--- allocate memory and initialize various variables
k=1;
img=sqrt(-1);
x_old=zeros(len1,1);
Nframes=floor(length(x)/len2)-1;
xfinal=zeros(Nframes*len2,1);
fr=Srate/nFFT;
fIndx=1:floor(4000/fr);
freqs=0:fr:(4000-fr);
%=============================== Start Processing =======================================================
%
k=1;
aa=0.98;
c=sqrt(pi)/2;
ksi_ar=zeros(Nframes,1);
gam_ar=ksi_ar;
C2=gamma(0.5);
%p=-1.0;
CC=gamma(0.6)/gamma(0.1);
STSA=0;
for n=1:Nframes
insign=win.*x(k:k+len-1);
%--- Take fourier transform of frame
spec=fft(insign,nFFT)/U;
sig=abs(spec); % compute the magnitude
sig2=sig.^2;
% numer = sqrt(gamma(1.5)*vk.*confhyperg(-0.5,1,-vk,40));
% denom = gammak.*sqrt(gamma(0.5)*confhyperg(0.5,1,-vk,40));
% hw=numer./denom;
k=k+len2;
end
%========================================================================================
if max(abs(xfinal))>32768
xfinal=xfinal*24000/max(abs(xfinal));
fprintf('Max amplitude exceeded 32768 for file %s\n',filename);
end
xfinal=xfinal/32768;
wavwrite(xfinal,Srate,16,outfile);
return;
% -------- plot SNR ---------------
%
timex=0:1000/Srate:Nframes*10;
subplot(2,1,1),h=plot(timex,x(1:length(timex)));
set(gca,'Box','off','Xlim',[0 timex(length(timex))],'FontSize',12);
ylabel('Amplitude');
time2=0:10:(Nframes-1)10;
subplot(2,1,2),
plot(time2,10log10(hgain));
h=plot(time2,10log10(ksi_ar),'b',time2,10log10(gam_ar),'r:');
legend('\xi_k','\gamma_k-1');
set(h,'Linewidth',1.5);
set(gca,'FontSize',12,'Box','off','Xlim',[0 time2(length(time2))],'Ylim',[-25 max(10*log10(gam_ar))]);
xlabel('Time (ms)');
ylabel(' SNR (dB)');
return;
subplot(3,1,3), %h=plot(time2,20*log10(hgain));
h=plot(20*log10(hgain));
%set(gca,'Box','off','Xlim',[0 timex(length(timex))],'FontSize',12);
ylabel('Amplitude');
function stsa_weuclid1(filename,outfile,p)
if nargin<3
fprintf('Usage: stsa_weuclid inFile outFile.wav p \n\n');
return;
end
fp=fopen(filename,'r');
if fp<=0, error('Unable to open input file'); end;
% ------------ get the file extension, and then read the whole file -----------
%
ind1=find(filename == '.');
if length(ind1)>1, ind=ind1(length(ind1)); else, ind=ind1; end;
ext = lower(filename(ind+1:length(filename)));
[HDRSIZE, Srate, bpsa, ftype] = gethdr(fp,ext);
x=fread(fp,inf,'short');
fclose(fp);
%x=x/32768;
% =============== Initialize variables ===============
%
len=floor(20Srate/1000); % Frame size in samples
if rem(len,2)==1, len=len+1; end;
PERC=50; % window overlap in percent of frame size
len1=floor(lenPERC/100);
len2=len-len1;
win=hanning(len); %tukey(len,PERC); % define window
U=norm(win);
% Noise magnitude calculations - assuming that the first 6 frames is noise/silence
%
nFFT=len;
nFFT2=len/2;
noise_mean=zeros(nFFT,1);
j=1;
for k=1:5
noise_mean=noise_mean+abs(fft(win.*x(j:j+len-1),nFFT)/U);
j=j+len;
end
noise_mu=noise_mean/5;
noise_mu2=noise_mu.^2;
%--- allocate memory and initialize various variables
k=1;
img=sqrt(-1);
x_old=zeros(len1,1);
Nframes=floor(length(x)/len2)-1;
xfinal=zeros(Nframes*len2,1);
fr=Srate/nFFT;
fIndx=1:floor(4000/fr);
freqs=0:fr:(4000-fr);
%=============================== Start Processing =======================================================
%
k=1;
aa=0.98;
c=sqrt(pi)/2;
ksi_ar=zeros(Nframes,1);
gam_ar=ksi_ar;
C2=gamma(0.5);
%p=-1.0;
CC=gamma(0.6)/gamma(0.1);
STSA=0;
for n=1:Nframes
insign=win.*x(k:k+len-1);
%--- Take fourier transform of frame
spec=fft(insign,nFFT)/U;
sig=abs(spec); % compute the magnitude
sig2=sig.^2;
gammak=min(sig2./noise_mu2,40); % post SNR if n==1 ksi=aa+(1-aa)*max(gammak-1,0); else ksi=aa*Xk_prev./noise_mu2 + (1-aa)*max(gammak-1,0); % a priori SNR %ksi=max(gammak-1,0); end vv=0.6; vk=ksi.*gammak./(vv+ksi); %----- weighted Euclidean distance ------------------------ numer=CC*sqrt(vk).*confhyperg(0.6,1,vk,100); denom=gammak.*confhyperg(0.1,1,vk,100); hw=numer./denom; %hw2=CC*sqrt(vk)./(gammak.*exp(-vk/2).*besseli(0,vk/2)); % if p=-1 % % ---- for the cosh measure
% numer = sqrt(gamma(1.5)*vk.*confhyperg(-0.5,1,-vk,40));
% denom = gammak.*sqrt(gamma(0.5)*confhyperg(0.5,1,-vk,40));
% hw=numer./denom;
% --- for the weighted cosh measure %CC2=gamma((p+3)/2)/gamma((p+1)/2); %numer=CC2*sqrt(vk).*sqrt(confhyperg(-(p+1)/2,1,-vk,40)); %denom=gammak.*sqrt(confhyperg(-(p-1)/2,1,-vk,40)); %hw=numer./denom; sig=sig.*hw; Xk_prev=sig.^2; xi_w= ifft( hw .* spec); xi_w= real( xi_w)*U; % --- Overlap and add --------------- % xfinal(k:k+ len2-1)= x_old+ xi_w(1:len1); x_old= xi_w(len1+ 1: len);
k=k+len2;
end
%========================================================================================
if max(abs(xfinal))>32768
xfinal=xfinal*24000/max(abs(xfinal));
fprintf('Max amplitude exceeded 32768 for file %s\n',filename);
end
xfinal=xfinal/32768;
wavwrite(xfinal,Srate,16,outfile);
return;
% -------- plot SNR ---------------
%
timex=0:1000/Srate:Nframes*10;
subplot(2,1,1),h=plot(timex,x(1:length(timex)));
set(gca,'Box','off','Xlim',[0 timex(length(timex))],'FontSize',12);
ylabel('Amplitude');
time2=0:10:(Nframes-1)10;
subplot(2,1,2),
plot(time2,10log10(hgain));
h=plot(time2,10log10(ksi_ar),'b',time2,10log10(gam_ar),'r:');
legend('\xi_k','\gamma_k-1');
set(h,'Linewidth',1.5);
set(gca,'FontSize',12,'Box','off','Xlim',[0 time2(length(time2))],'Ylim',[-25 max(10*log10(gam_ar))]);
xlabel('Time (ms)');
ylabel(' SNR (dB)');
return;
subplot(3,1,3), %h=plot(time2,20*log10(hgain));
h=plot(20*log10(hgain));
%set(gca,'Box','off','Xlim',[0 timex(length(timex))],'FontSize',12);
ylabel('Amplitude');
相关文章推荐
- ubuntu 安装 matlab
- TeX系列: MATLAB和LaTeX结合绘图
- MATLAB中的bwlookup函数
- matlab窗口输出函数
- MATLAB 基础之程序设计
- Matlab xml读写
- matlab 修改xml文件的属性
- A release of a Camera Calibration Toolbox for Matlab
- Camera Calibration Toolbox for Matlab
- Matlab2013 配置vlfeat
- 3D散点 表面,2D散点,边界。 matlab 儿子的papa
- 凸优化工具箱cvx
- 主成分分析法原理与MATLAB实现
- matlab—多行注释
- MATLAB 人机交互(窗口)设置命令 (转载)
- Matlab中plot函数全功能解析
- matlab xlswrite 用法以及实例(转自“百度经验”)
- matlab中导入excel中的数据(转自“百度经验”)——亲测可用
- matlab 实现图像PSNR的小程序
- Matlab中的lsqcurvefit,非线性拟合