用matlab编程设计一个巴特沃斯滤波器，对声音进行去噪

fs=44100

[x,fs,bits]=wavread('ding11.wav')

%sound(x)

t=0:(size(x)-1)

x2=rand(1,length(x))' %产生一与x长高裤州度一致的随机信号

y=x+x2

%加入正弦噪声

t=0:(n-1)

Au=0.03

d=[Au*sin(2*pi*500*t)]'

y=x+d

wp=0.25*pi

ws=0.3*pi

wdelta=ws-wp

N=ceil(6.6*pi/wdelta) %取整

wn=(0.2+0.3)*pi/2

b=fir1(N,wn/pi,hamming(N+1)) %选择窗函数，并归一化截止频率

figure(1)

freqz(b,1,512)

f2=filter(bz,az,y)

figure(2)

subplot(2,1,1)

plot(t,y)

title('滤波前的时域波形')

subplot(2,1,2)

plot(t,f2)

title('滤波后的时域波形戚蔽')

sound(f2) %播放滤波后的语音信号

F0=fft(f1,1024)

f=fs*(0:511)/1024

figure(3)

y2=fft(y,1024)

subplot(2,1,1)

plot(f,abs(y2(1:512)))%画出滤波前的频谱图

title('滤波前纯坦的频谱')

xlabel('Hz')

ylabel('fuzhi')

subplot(2,1,2)

F1=plot(f,abs(F0(1:512))) %画出滤波后的频谱图

title('滤波后的频谱')

xlabel('Hz')

ylabel('fuzhi')

fs = 44100 %采样率

f0 = 5000   %信号频率

N = 1024

%巴特沃斯低通滤波器

Wp = 10000/fs 

Ws = 15000/fs

Rp = 3

Rs = 60

[n,Wn] = buttord(Wp,Ws,Rp,Rs)

[b,a] = butter(n,Wn)

figure

freqz(b,a,N) 

title('巴特沃斯汪铅派低通滤波器特性')

tp = N/fs %采样时长困贺

t = 0:1/fs:tp

y = sin(2*pi*f0*t)  %信号

yn = y + rand(1,N+1) %加噪声

%显示10个周期

t2 = 0:1/fs:10/f0

L = length(t2)

figure

subplot(311)plot(t2,y(1:L))title('信号')ylim([-2,2])

subplot(312)plot(t2,yn(1:L))title('信号加噪声激弯')ylim([-2,2])

%滤波

yf = filter(b,a,yn)

subplot(313)plot(t2,yf(1:L))title('滤波后信号')ylim([-2,2])

winsize=256%窗长

n=0.1%噪声水平

a=2

b=6

[speech,fs,nbits]=wavread('E:\matlab\louyin.wav')%读入wav文件

size=length(speech)%语音长度

numofwin=floor(size/winsize)%窗举斗数

%定义汉明窗

ham=hamming(winsize)'

hamwin=zeros(1,size)

enhanced=zeros(1,size)

improved=zeros(1,size)

%生成噪声信号

noise=n*randn(1,size)

y=speech'+noise

%噪声处理

noisy=n*randn(1,winsize)

N=fft(noisy)

npow=abs(N)

for q=1:2*numofwin-1

yframe=y(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)%分正正磨帧

hamwin(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)=hamwin(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)+ham%

%加噪信号FFT

y1=fft(yframe.*ham)

ypow=abs(y1)%加噪信号幅度

yangle=angle(y1)%相位

%计算功率谱密度

Py=ypow.^2

Pn=npow.^2

Pyy=ypow.^a

Pnn=npow.^a

%基本谱减

for i=1:winsize

if Py(i)-Pn(i)>0

Ps(i)=Py(i)-Pn(i)

else

Ps(i)=0

end

end

s=sqrt(Ps).*exp(j*yangle)

for i=1:winsize

if Pyy(i)-b*Pnn(i)>0

Pss(i)=Pyy(i)-b*Pnn(i)

else

Pss(i)=0

end

end

ss=Pss.^(1/a).*exp(j*yangle)

%去噪语音IFFT

enhanced(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)=enhanced(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)+real(ifft(s))

improved(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)=improved(1+(q-1)*winsize/2:winsize+(q-1)*winsize/2)+real(ifft(ss))

end

%去除汉明窗引起的增益

for i=1:size

if hamwin(i)==0

enhanced(i)=0

improved(i)=0

else

enhanced(i)=enhanced(i)/hamwin(i)

improved(i)=improved(i)/hamwin(i)

end

end

SNR1=10*log10(var(speech')/清或var(noisy))%加噪语音信噪比

SNR2=10*log10(var(speech')/var(enhanced-speech'))%增强语音信噪比

SNR3=10*log10(var(speech')/var(improved-speech'))

figure(1)plot(speech')%原始语音波形

title(['Original Voice(n=',num2str(n),')'])

figure(2)plot(y)

title(['Noise Added(SNR=',num2str(SNR1),'dB)'])

figure(3)plot(enhanced)

title(['Enhanced Voice(SNR=',num2str(SNR2),'dB)'])

figure(4)plot(improved)

title(['Improved Voice(SNR=',num2str(SNR3),'dB)'])

---