语音信号PCM编译码 MATLAB程序

% clc

% clear all

% 产生信号

load handel

x=y(1:20000)%取前20000个采样点

sound(x,Fs)

%PCM编码

x1=x/0.8.*2048

yy=pcm_encode(x1)

figure

subplot(2,1,1)

stem(yy(1:80),'.')

title('PCM编码后的波形')

%加噪声

snr=10

sp=mean(yy.^2)

attn=sp./ 10^(snr/10)

attn = sqrt(attn)

noise=randn(1,length(yy)).*attn

np=mean(noise.^2)

snr1=10*log10(sp/np)

data=yy+noise

% data=yy%不加噪声

subplot(2,1,2)

stem(data(1:80),'羡厅和.r')

title('PCM加噪声后波形')

%译码

demodata=data>0.5

zz=pcm_decode(demodata,0.8)

figure

subplot(2,1,1)

plot(x)

title('兄盯原始语音信号')

subplot(2,1,2)

plot(zz)

title('译码后的语音信号')

sound(zz,Fs)

figure

plot(x,'b')

hold on

plot(zz,'r')

legend('编译码前的语音信号','编译码后的语音信号')

title('编，译码前后的语音信号')

function y = pcm_encode( x )

y=zeros(length(x),8) %存储矩阵（全零）

z=sign(x) %判断x的正负

x=abs(x)%取绝伏仿对值

%%段落码判断段区间的取值范围为前开后闭区间

for k=1:length(x)

%符号位的判断

if z(k)>0

y(k,1)=1

elseif z(k)<0

y(k,1)=0

end

if x(k)>128 &x(k)<=2048 %在第五段与第八段之间，段位码第一位都为“1”

y(k,2)=1

end

if (x(k)>32 &x(k)<=128) || (x(k)>512 &x(k)<=2048)

y(k,3)=1 %在第三四七八段内，段位码第二位为“1”

end

if (x(k)>16&x(k)<=32)||(x(k)>64&x(k)<=128)||(x(k)>256&x(k)<=512)||(x(k)>1024&x(k)<=2048)

y(k,4)=1 %在二四六八段内，段位码第三位为“1”

end

end

%段内码判断程序

N=zeros(1,length(x))

for k=1:length(x)

N(k)=y(k,2)*4+y(k,3)*2+y(k,4)+1 %找到x位于第几段

end

a=[0,16,32,64,128,256,512,1024] %量化间隔

b=[1,1,2,4,8,16,32,64]%除以16，得到每段的最小量化间隔

for m=1:length(x)

q=ceil((x(m)-a(N(m)))/b(N(m))) %求出在段内的位置

if q==0

y(m,(5:8))=[0,0,0,0] %如果输入为零则输出“0”

else k=num2str(dec2bin(q-1,4)) %编码段内码为二进制

y(m,5)=str2num(k(1))

y(m,6)=str2num(k(2))

y(m,7)=str2num(k(3))

y(m,8)=str2num(k(4))

end

end

%将N行8列矩阵转换为1行8*N列的矩阵

y=y'

y=reshape(y,1,length(x)*8)

end

function x=pcm_decode(y,max)

%将1行8*N列的矩阵转换为N行8列矩阵

y=reshape(y,8,length(y)/8)

y=y'

%PCM译码

n=size(y,1) %求出输入码组的个数

a=[0,16,32,64,128,256,512,1024] %段落起点值

b=[1,1,2,4,8,16,32,64] %每段的最小量化间隔

for k=1:n

t1=y(k,1)%取符号

t2=y(k,2)*4+y(k,3)*2+y(k,4)+1%判断段落位置

t3=y(k,5)*8+y(k,6)*4+y(k,7)*2+y(k,8) %判断段内位置

if t3==0 %段内码为零时

m(k)=(a(t2)+1+0.5*b(t2))/2048*max

else

m(k)=(a(t2)+b(t2)*t3+0.5*b(t2))/2048*max%还原出量化后的电平值

end

%判断符号位

if t1==0

x(k)=-m(k)

else

x(k)=m(k)

end

end

end

不好意思，没看到维纳滤波，程序就不删了

%谱减法语音增强

%输入参数s 语音数据，fs 采样频率，p 下面有说明，共11个，可不差正输入，有默认值

%“过度减法(oversubtraction)”作减法的时候,保留一小部分原来的背景噪音,用这部分背景噪音来掩盖住音乐噪音虚行悔的谱峰,从而消除了令人不悦的音乐噪音。

%通过给的参数p，估计噪音，做谱减法。从而消除噪音。

function [ss,po]=specsubm(s,fs,p)

%利用频谱相减（spectral subtraction）增强 [SS,PO]=(S,FS,P)

%

% implementation of spectral subtraction algorithm by R Martin (rather slow)

% algorithm parameters: t* in seconds, f* in Hz, k* dimensionless

% 1: tg = smoothing time constant for signal power estimate (0.04): high=reverberant, low=musical

% 2: ta = smoothing time constant for signal power estimate

%used in noise estimation (0.1)

% 3: tw = fft window length (will be rounded up to 2^nw samples)

% 4: tm = length of minimum filter (1.5): high=slow response to noise increase, low=distortion

% 5: to = time constant for oversubtraction factor (0.08)

% 6: fo = oversubtraction corner frequency (800): high=distortion, low=musical

% 7: km = number of minimisation buffers to use (4): high=waste memory, low=noise modulation

% 8: ks = oversampling constant (4)

% 9: kn = noise estimate compensation (1.5)

% 10:kf = subtraction floor (0.02): high=noisy, low=musical

% 11:ko = oversubtraction scale factor (4): high=distortion, low=musical

%检查函数带森的输入参数，如果输入少于三个，po为默认值，po的参数上面有说明

if nargin<3 po=[0.04 0.1 0.032 1.5 0.08 400 4 4 1.5 0.02 4].'else po=pend

ns=length(s)

ts=1/fs

ss=zeros(ns,1)

ni=pow2(nextpow2(fs*po(3)/po(8)))

ti=ni/fs

nw=ni*po(8)

nf=1+floor((ns-nw)/ni)

nm=ceil(fs*po(4)/(ni*po(7)))

win=0.5*hamming(nw+1)/1.08win(end)=[]

zg=exp(-ti/po(1))

za=exp(-ti/po(2))

zo=exp(-ti/po(5))

px=zeros(1+nw/2,1)

pxn=px

os=px

mb=ones(1+nw/2,po(7))*nw/2

im=0

osf=po(11)*(1+(0:nw/2).'*fs/(nw*po(6))).^(-1)

imidx=[13 21]'

x2im=zeros(length(imidx),nf)

osim=x2im

pnim=x2im

pxnim=x2im

qim=x2im

for is=1:nf

idx=(1:nw)+(is-1)*ni

x=rfft(s(idx).*win)

x2=x.*conj(x)

pxn=za*pxn+(1-za)*x2

im=rem(im+1,nm)

if im

mb(:,1)=min(mb(:,1),pxn)

else

mb=[pxn,mb(:,1:po(7)-1)]

end

pn=po(9)*min(mb,[],2)

%os= oversubtraction factor

os=zo*os+(1-zo)*(1+osf.*pn./(pn+pxn))

px=zg*px+(1-zg)*x2

q=max(po(10)*sqrt(pn./x2),1-sqrt(os.*pn./px))

ss(idx)=ss(idx)+irfft(x.*q)

end

if nargout==0

soundsc([sss],fs)

end

---