摘要
一、引言
二、正文
1。設計要求
2。設計步驟
3。設計內容
4。簡易GUI設計
三、結論
四、收穫與心得
五、附錄
一、引言
隨著Matlab模擬技術的推廣，我們可以在計算機上對聲音訊號進行處理，甚至是模擬。透過計算機作圖，取樣，我們可以更加直觀的瞭解語音訊號的性質，透過matlab程式設計，呼叫相關的函式，我們可以非常方便的對訊號進行運算和處理。
二、正文
2。1 設計要求
在有噪音的環境中錄製語音，並設計濾波器去除噪聲。
2。2 設計步驟
分析原始訊號，畫出原始訊號頻譜圖及時頻圖，確定濾波器型別及相關指標 ；
按照型別及指標要求設計出濾波器，畫出濾波器幅度和相位響應，分析該濾波器是否符合要求；
用所設計的濾波器對原始訊號進行濾波處理，畫出濾波後訊號的頻譜圖及時頻圖；
對濾波前的訊號進行分析比對，評估所設計濾波器效能。
2。3 設計內容
1。原始訊號分析
分析訊號的譜圖可知，噪音在1650HZ和3300HZ附近的能量較高，而人聲的能量基本位於1000HZ以下。因此，可以設計低通濾波器對訊號進行去噪處理。
2。IIR濾波器設計
用雙線性變換法分別設計了巴特沃斯低通濾波器和橢圓低通濾波器和帶阻濾波器：
①巴特沃斯濾波器
fp=800；fs=1300；rs=35；rp=0。5；
程式程式碼如下：
fp=800；fs=1300；rs=35；rp=0。5；Fs=44100；
wp=2*Fs*tan（2*pi*fp/（2*Fs））；ws=2*Fs*tan（2*pi*fs/（2*Fs））；
［n，wn］=buttord（wp，ws，rp，rs，‘s’）；
［b，a］=butter（n，wn，‘s’）；
［num，den］=bilinear（b，a，Fs）；
［h，w］=freqz（num，den，512，Fs）；
②橢圓低通濾波器
fp=1300；fs=1600；rs=60；rp=0。5；
程式程式碼如下：
fp=1300；fs=1600；rs=60；rp=0。5；Fs=44100；
wp=2*Fs*tan（2*pi*fp/（2*Fs））；ws=2*Fs*tan（2*pi*fs/（2*Fs））；
［n，wn］=ellipord（wp，ws，rp，rs，‘s’）；
［b，a］=ellip（n，rp，rs，wn，‘s’）；
［num，den］=bilinear（b，a，Fs）；
［h，w］=freqz（num，den，512，Fs）；
③帶阻濾波器
fp1=800；fp2=2300；fs1=1300；fs2=1800；rs=30；rp=0。6
fp3=2800；fp4=4000；fs3=3200；fs4=3700；rs=30；rp=0。6
程式程式碼如下：
fp1=800；fp2=2300；fs1=1300；fs2=1800；rs=30；rp=0。6；Fs=44100；
fp=［fp1，fp2］；fs=［fs1，fs2］；
wp=2*Fs*tan（2*pi*fp/（2*Fs））；
ws=2*Fs*tan（2*pi*fs/（2*Fs））；
［n，wn］=buttord（wp，ws，rp，rs，‘s’）；
［b，a］=butter（n，wn，‘stop’，‘s’）；
［num，den］=bilinear（b，a，Fs）；
［h，w］=freqz（num，den，512，Fs）；
fp3=2800；fp4=4000；fs3=3200；fs4=3700；rs=30；rp=0。6；Fs=44100；
fp1=［fp3，fp4］；fs1=［fs3，fs4］；
wp1=2*Fs*tan（2*pi*fp1/（2*Fs））；
ws1=2*Fs*tan（2*pi*fs1/（2*Fs））；
［n1，wn1］=buttord（wp1，ws1，rp，rs，‘s’）；
［b1，a1］=butter（n1，wn1，‘stop’，‘s’）；
［num1，den1］=bilinear（b1，a1，Fs）；
［h1，w1］=freqz（num1，den1，512，Fs）；
FIR濾波器
①加hamming窗
n=100；fp=1000；Fs=44100；
b=fir1（n，fp/（Fs/2），Hamming（n+1））；
［h，w］=freqz（b，1，512，Fs）；
②加hanning窗
n=；fp=1000；Fs=44100；
b=fir1（n，fp/（Fs/2），Hanning（n+1））；
［h，w］=freqz（b，1，512，Fs）；
③加blackman窗
n=100；fp=1000；Fs=44100；
b=fir1（n，fp/（Fs/2），blackman（n+1））；
［h，w］=freqz（b，1，512，Fs）；
濾波前後比對
①巴特沃斯低通濾波器濾波後
②橢圓低通濾波器濾波後
③帶阻濾波器
④加hamming窗
⑤加hanning窗
⑥加blackman窗
2。4簡易GUI介面設計
為了便於操作和演示，設計瞭如下的簡易GUI介面。
結論
由以上譜圖分析可知，經過濾波器濾波後，訊號中的高頻雜音明顯被抑制，而人聲成分大部分被保留，起到了預期的濾波作用。
對比所設計的兩種濾波器，橢圓濾波器在過渡帶相對較窄的情況下，能滿足相對較高阻帶衰減。
四、收穫與心得
本次設計大概進行了五週的時間，語音訊號處理的是目前比較流行且十分有趣的，在程式設計實現的過程中還是遇到了很多困難。在前期的準備工作中，查閱了大量資料，以完善我們的理論知識。我們為了完成本次設計，我們透過查閱相關書籍以及matlab中的幫助，選用不同的matlab函式，嘗試不同的引數。經過接近一個禮拜的反覆除錯，最終基本的實現了設計任務。
雖然遇到了很多困難，但是我們在設計過程中都有收穫很大。本次設計將訊號與系統課上學習的知識用於實踐，讓我們對對語音訊號處理更深入的瞭解，也讓我們加深了對濾波器相關內容的理解，同時也使得我們的Matlab能力有了很大的提高。
參考文獻
《應用matlab實現訊號分析和處理》 科學出版社
附錄
1 巴特沃斯低通濾波器
fp=800；fs=1300；rs=35；rp=0。5；Fs=44100；
wp=2*Fs*tan（2*pi*fp/（2*Fs））；
ws=2*Fs*tan（2*pi*fs/（2*Fs））；
［n，wn］=buttord（wp，ws，rp，rs，‘s’）；
［b，a］=butter（n，wn，‘s’）；
［num，den］=bilinear（b，a，Fs）；
［h，w］=freqz（num，den，512，Fs）；
figure（1）
%subplot（3，1，1）
plot（w，abs（h））；
xlabel（‘頻率/Hz’）；ylabel（‘幅值’）；
title（‘巴特沃斯低通濾波器幅度特性’）；
axis（［0，5000，0，1。2］）
grid on；
figure（2）
%subplot（3，1，2）
plot（w，20*log10（abs（h）））；
xlabel（‘頻率/Hz’）；ylabel（‘幅值db’）；
title（‘巴特沃斯低通濾波器幅度特性db’）；
axis（［0，5000，-90，10］）；
grid on；
figure（3）
plot（w，180/pi*unwrap（angle（h）））；
xlabel（‘頻率/Hz’）；ylabel（‘相位’）；
title（‘巴特沃斯低通濾波器相位特性’）；
axis（［0，5000，-1000，10］）
grid on；
［s1，Fs，bits］=wavread（‘D：\222。wav’）；
x1=s1（：，1）；sound（x1，Fs，bits）；
N1=length（x1）；
Y1=fft（x1，N1）；
f1=Fs*（0：N1-1）/N1；
t1=（0：N1-1）/Fs；
figure（4）
plot（f1，abs（Y1））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘原始訊號頻譜’）；
grid on；axis（［0 6000 0 400］）
y=filter（num，den，x1）；
sound（y，Fs，bits）；
N2=length（y）；
Y2=fft（y，N2）；
f2=Fs*（0：N2-1）/N2；
t2=（0：N2-1）/Fs；
figure（5）
plot（f2，abs（Y2））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘過濾後訊號的頻譜’）；
grid on；axis（［0 6000 0 100］）
2橢圓低通濾波器
fp=1300；fs=1600；rs=60；rp=0。5；Fs=44100；
wp=2*Fs*tan（2*pi*fp/（2*Fs））；
ws=2*Fs*tan（2*pi*fs/（2*Fs））；
［n，wn］=ellipord（wp，ws，rp，rs，‘s’）；
［b，a］=ellip（n，rp，rs，wn，‘s’）；
［num，den］=bilinear（b，a，Fs）；
［h，w］=freqz（num，den，512，Fs）；
figure（1）
plot（w，abs（h））；
xlabel（‘頻率/Hz’）；ylabel（‘幅值’）；
title（‘橢圓低通濾波器幅度特性’）；
axis（［0，5000，0，1。2］）
grid on；
figure（2）
plot（w，20*log10（abs（h）））；
xlabel（‘頻率/Hz’）；ylabel（‘幅值db’）；
title（‘橢圓低通濾波器幅度特性db’）；
axis（［0，5000，-90，10］）；
grid on；
figure（3）
plot（w，180/pi*unwrap（angle（h）））；
xlabel（‘頻率/Hz’）；ylabel（‘相位’）；
title（‘橢圓低通濾波器相位特性’）；
axis（［0，5000，-1000，10］）
grid on；
［s1，Fs，bits］=wavread（‘D：\222。wav’）；
x1=s1（：，1）；sound（x1，Fs，bits）；
N1=length（x1）；
Y1=fft（x1，N1）； %對訊號做N點FFT變換
f1=Fs*（0：N1-1）/N1；
t1=（0：N1-1）/Fs；
figure（4）
plot（f1，abs（Y1））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘原始訊號頻譜’）；
grid on；axis（［0 6000 0 400］）
y=filter（num，den，x1）；
sound（y，Fs，bits）；
N2=length（y）；
Y2=fft（y，N2）； %對訊號做N點FFT變換
f2=Fs*（0：N2-1）/N2；
t2=（0：N2-1）/Fs；
figure（5）
plot（f2，abs（Y2））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘過濾後訊號的頻譜’）；
grid on；axis（［0 6000 0 100］）
3.帶阻濾波器
fp1=800；fp2=2300；fs1=1300；fs2=1800；rs=30；rp=0。6；Fs=44100；
fp=［fp1，fp2］；fs=［fs1，fs2］；
wp=2*Fs*tan（2*pi*fp/（2*Fs））；
ws=2*Fs*tan（2*pi*fs/（2*Fs））；%wap=2*tan（wp/2）/Ts
［n，wn］=buttord（wp，ws，rp，rs，‘s’）；
［b，a］=butter（n，wn，‘stop’，‘s’）；
［num，den］=bilinear（b，a，Fs）；
［h，w］=freqz（num，den，512，Fs）；
fp3=2800；fp4=4000；fs3=3200；fs4=3700；rs=30；rp=0。6；Fs=44100；
fp1=［fp3，fp4］；fs1=［fs3，fs4］；
wp1=2*Fs*tan（2*pi*fp1/（2*Fs））；
ws1=2*Fs*tan（2*pi*fs1/（2*Fs））；%wap=2*tan（wp/2）/Ts
［n1，wn1］=buttord（wp1，ws1，rp，rs，‘s’）；
［b1，a1］=butter（n1，wn1，‘stop’，‘s’）；
［num1，den1］=bilinear（b1，a1，Fs）；
［h1，w1］=freqz（num1，den1，512，Fs）；
figure（1）
plot（w，abs（h））；
xlabel（‘頻率/Hz’）；ylabel（‘幅值’）；
title（‘巴特沃斯帶阻濾波器幅度特性’）；
axis（［0，5000，0，1。2］）
grid on；
figure（2）
plot（w，20*log10（abs（h）））；
xlabel（‘頻率/Hz’）；ylabel（‘幅值db’）；
title（‘巴特沃斯帶阻濾波器幅度特性db’）；
axis（［0，5000，-90，10］）；
grid on；
figure（3）
plot（w，180/pi*unwrap（angle（h）））；
xlabel（‘頻率/Hz’）；ylabel（‘相位’）；
title（‘巴特沃斯帶阻濾波器相位特性’）；
axis（［0，5000，-1000，10］）
grid on；
figure（4）
plot（w1，abs（h1））；
xlabel（‘頻率/Hz’）；ylabel（‘幅值’）；
title（‘巴特沃斯帶阻濾波器幅度特性’）；
axis（［0，5000，0，1。2］）
grid on；
figure（5）
plot（w1，20*log10（abs（h1）））；
xlabel（‘頻率/Hz’）；ylabel（‘幅值db’）；
title（‘巴特沃斯帶阻濾波器幅度特性db’）；
axis（［0，5000，-90，10］）；
grid on；
figure（6）
plot（w1，180/pi*unwrap（angle（h1）））；
xlabel（‘頻率/Hz’）；ylabel（‘相位’）；
title（‘巴特沃斯帶阻濾波器相位特性’）；
axis（［0，5000，-1000，10］）
grid on；
［s1，Fs，bits］=wavread（‘D：\222。wav’）；
x1=s1（：，1）；sound（x1，Fs，bits）；
N1=length（x1）；
Y1=fft（x1，N1）； %對訊號做N點FFT變換
f1=Fs*（0：N1-1）/N1；
t1=（0：N1-1）/Fs；
figure（7）
plot（f1，abs（Y1））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘原始訊號頻譜’）；
grid on；axis（［0 6000 0 400］）
y1=filter（num，den，x1）；
y=filter（num1，den1，y1）；
sound（y，Fs，bits）；
N2=length（y）；
Y2=fft（y，N2）；
f2=Fs*（0：N2-1）/N2；
t2=（0：N2-1）/Fs；
figure（8）
plot（f2，abs（Y2））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘過濾後訊號的頻譜’）；
grid on；axis（［0 6000 0 100］）
加hamming窗
n=100；fp=1000；Fs=44100；
b=fir1（n，fp/（Fs/2），Hamming（n+1））；
［h，w］=freqz（b，1，512，Fs）；
a=num2str（a）；
［s1，Fs，bits］=wavread（‘D：\222。wav’）；
x1=s1（：，1）；
sound（x1，Fs，bits）；
N1=length（x1）；
Y1=fft（x1，N1）；
f1=Fs*（0：N1-1）/N1；
y=fftfilt（b，x1）；
sound（y，Fs，bits）；
Y2=fft（y，N2）；
f2=Fs*（0：N2-1）/N2；
figure（4）
subplot（2，1，1）
plot（f1，abs（Y1））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘原始訊號頻譜’）；
grid on；axis（［0 6000 0 400］）
subplot（2，1，2）
plot（f2，abs（Y2））；
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘過濾後訊號的頻譜’）；
grid on；axis（［0 6000 0 100］）
加hanning窗
n=100；fp=1000；Fs=44100；
b=fir1（n，fp/（Fs/2），Hanning（n+1））；
［h，w］=freqz（b，1，512，Fs）；
a=num2str（a）；
［s1，Fs，bits］=wavread（‘D：\222。wav’）；
x1=s1（：，1）；
N1=length（x1）；
Y1=fft（x1，N1）；
f1=Fs*（0：N1-1）/N1；
y=fftfilt（b，x1）；
sound（y，Fs，bits）；
N2=length（y）；
Y2=fft（y，N2）；
f2=Fs*（0：N2-1）/N2；
figure（4）
subplot（2，1，1）
plot（f1，abs（Y1））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘原始訊號頻譜’）；
grid on；axis（［0 6000 0 400］）
subplot（2，1，2）
plot（f2，abs（Y2））；
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘過濾後訊號的頻譜’）；
grid on；axis（［0 6000 0 100］）
加blackman窗
n=100；fp=1000；Fs=44100；
b=fir1（n，fp/（Fs/2），blackman（n+1））；
［h，w］=freqz（b，1，512，Fs）；
a=num2str（a）；
［s1，Fs，bits］=wavread（‘D：\222。wav’）；
x1=s1（：，1）；
sound（x1，Fs，bits）；
N1=length（x1）；
Y1=fft（x1，N1）；
f1=Fs*（0：N1-1）/N1；
y=fftfilt（b，x1）；
sound（y，Fs，bits）；
N2=length（y）；
f2=Fs*（0：N2-1）/N2；
figure（4）
subplot（2，1，1）
plot（f1，abs（Y1））
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘原始訊號頻譜’）；
grid on；axis（［0 6000 0 400］）
subplot（2，1，2）
plot（f2，abs（Y2））；
xlabel（‘頻率/Hz’）；ylabel（‘幅度’）；
title（‘過濾後訊號的頻譜’）；
grid on；axis（［0 6000 0 100］）
7。巴特沃斯低通濾波時頻分析
clear
Fs=44100；
［s1，Fs，bits］=wavread（‘D：\222。wav’）；
x=s1（：，1）；x1=x（40001：64576）
l=length（x1）；
N=6；
m=downsample（x1，N）；%降抽樣後 Fs=7350hz
z=length（m）；%4096
［tfr， t， f］ = tfrstft（m）；
figure（13）
contour（t，（Fs/N）*（0：z/2-1）/z，abs（tfr（1：z/2，：））。^2）；
xlabel（‘時間t’）；
ylabel（‘頻率f’）；
title（‘等高線圖（濾波前）’）；
grid on；
figure（14）
mesh（t，（Fs/N）*（0：z/2-1）/z，abs（tfr（1：z/2，：））。^2）；
xlabel（‘時間t’）；
ylabel（‘頻率f’）；
zlabel（‘幅值A’）；
title（‘三維圖（濾波前）’）；
grid on；
fp=1300；fs=1600；rs=60；rp=0。5；Fs=44100；
wp=2*Fs*tan（2*pi*fp/（2*Fs））；ws=2*Fs*tan（2*pi*fs/（2*Fs））；
［n，wn］=ellipord（wp，ws，rp，rs，‘s’）；
［b，a］=ellip（n，rp，rs，wn，‘s’）；
［num，den］=bilinear（b，a，Fs）；
［h，w］=freqz（num，den，512，Fs）；
y=filter（num，den，x）；
x1=y（40001：64576）
l=length（x1）；
N=6；
m=downsample（x1，N）；%降抽樣後 Fs=7350hz
z=length（m）；%4096
［tfr，t，f］ = tfrstft（m）；
figure（15）
contour（t，（Fs/N）*（0：z/2-1）/z，abs（tfr（1：z/2，：））。^2）；
xlabel（‘時間t’）；
ylabel（‘頻率f’）；
title（‘等高線圖（濾波後）’）；
grid on；
figure（16）
mesh（t，（Fs/N）*（0：z/2-1）/z，abs（tfr（1：z/2，：））。^2）；
xlabel（‘時間t’）；
ylabel（‘頻率f’）；
zlabel（‘幅值A’）；
title（‘三維圖（濾波後）’）；
grid on；
