數字信號Matlab作業(yè)

1、生物醫(yī)學信號處理大作業(yè)題目：布萊克曼窗設計FIR濾波器處理語音信號 學生姓名 學院名稱 精密儀器與光電子工程 專 業(yè) 生物醫(yī)學工程?班 學 號 3011202312 題目要求： 錄制自己的一段語音：“天津大學精密儀器與光電子工程學院生物醫(yī)學工程X班XXX, College of precision instrument and opto-electronics engineering, biomedical engineering”，時間控制在15秒到30秒左右；利用wavread函數對自己的語音進行采樣，記住采樣頻率。 （1）求原始語音信號的特征頻帶（比如譜峰位置）：可以分別對一定時間間隔內

2、，求功率譜（傅里葉變換結果取模的平方）并畫出功率譜。 （2）根據語音信號頻譜特點，設計FIR或IIR濾波器，分別畫出濾波器幅頻和相頻特性曲線。說明濾波器特性參數。用設計的濾波器對信號濾波，畫出濾波前后信號的頻譜圖。用sound函數回放語音信號,說明利用高通/低通/帶通濾波后的效果，不同特征頻帶被濾除后分別有什么效果。 （3）求出特征頻段語音信號隨時間變化的曲線（每隔一定時間求一次功率譜，連接成曲線，即短時 FFT）。 （4）選做：語譜圖：橫軸為時間，縱軸為頻率，灰度值大小表示功率譜值的大小。(提示，可以采用spectrogram函數) (5) 選做：分析自己的語音頻譜特點，比如中英文發(fā)音的區(qū)別

3、。1. 求原始語音信號的特征頻帶代碼：clcclear allFile=('I:Matlab大作業(yè)work.wav')Y1,Fs,nbits=wavread(File);Y1 = Y1(:,1)'%取雙聲道中的一個聲道l=length(Y1);Y=Y1(1:l);Fs=44100;figure(1),plot(Y);title('原聲音信號波形') for i = 1:18; x = Y(Fs*(i-1)+1 : Fs*i); z = fft(x , Fs); p = z .*conj(z)/Fs; %功率譜 figure(i); plot(p); xl

4、abel('頻率Hz') ylabel('功率') axis(0 2500 0 1);end1.1 原始語音信號1.2 以1秒為時間間隔，求出17秒的語音信號的功率譜如下：2. 根據語音信號頻譜特點，設計FIR或IIR濾波器，分別畫出濾波器幅頻和相頻特性曲線。說明濾波器特性參數。用設計的濾波器對信號濾波，畫出濾波前后信號的頻譜圖。用sound函數回放語音信號,說明利用高通/低通/帶通濾波后的效果，不同特征頻帶被濾除后分別有什么效果。代碼：%=低通濾波器=%fp=1300; fs=1500; % 輸入給定指標% 用窗函數法設計濾波器wc=(fp+fs)/Fs; %

5、理想低通濾波器截止頻率(關于pi歸一化）B=2*pi*(fs-fp)/Fs; %過渡帶寬度指標N=ceil(12*pi/B); %blackman窗的長度Nh1=fir1(N-1,wc,'low',blackman(N);%用窗函數法設計法的繪圖部分h11,w1=freqz(h1,1);figure(2)subplot(2,1,1);plot(w1*Fs/pi/2,20*log10(abs(h11);title('低通濾波器');grid on;axis(0,2500,-150,10);xlabel('頻率/Hz');ylabel('幅度

6、/dB')subplot(2,1,2)plot(w1*Fs/pi/2,angle(h11);grid on;axis(0,1000,-4,4);xlabel('頻率/Hz');ylabel('角度/rad');%sound(ywt,Fs); %=高通濾波器=%fp=300; fs=100; % 輸入給定指標%窗函數法設計濾波器wc=(fp+fs)/Fs; %理想高通濾波器截止頻率(關于pi歸一化）B=2*pi*(fp-fs)/Fs; %過渡帶寬度指標N=ceil(12*pi/B); %blackman窗的長度Nh2=fir1(N-1,wc,'hi

7、gh',blackman(N);%=以下為繪圖部分=h22,w2=freqz(h2,1);figure(3)subplot(2,1,1);plot(w2*Fs/pi/2,20*log10(abs(h22);title('高通濾波器');grid on;axis(0,1000,-100,20);xlabel('頻率/Hz');ylabel('幅度/dB')subplot(2,1,2)plot(w2*Fs/pi/2,angle(h22);grid on;axis(0,1000,-4,4);xlabel('頻率/Hz');yla

8、bel('角度/rad'); %=語音信號經過低通和高通濾波器=Y_L=filter(h1,1,Y);%sound(Y_L,Fs);Y_H=filter(h2,1,Y);%sound(Y_H,Fs);Y_N=filter(h2,1,Y_L);%sound(Y_N,Fs); %=原始語音信號的頻譜圖=DFT = fft(Y,Fs);DFT = DFT./max(DFT);figure(4),plot(abs(DFT);axis(0 10000 0 0.5);grid; %=語音信號經過低通濾波器后的頻譜圖=DFT_L = fft(Y_L,Fs);DFT_L = DFT_L./ma

9、x(DFT_L);figure(5),plot(abs(DFT_L);axis(0 2000 0 0.8);grid; %=語音信號經過高通濾波器后的頻譜圖=DFT_H = fft(Y_H,Fs);DFT_H = DFT_H./max(DFT_H);figure(6),plot(abs(DFT_H);axis(0 8000 0 0.7);grid; %=語音信號經過帶通濾波器后的頻譜圖=DFT_N = fft(Y_N,Fs);DFT_N = DFT_N./max(DFT_N);figure(7),plot(abs(DFT_N);axis(0 2000 0 1);grid; （1） 用布萊克曼窗

10、FIR濾波器，畫出低通濾波器的幅頻和相頻特性曲線如下： 低通濾波器的幅頻和相頻特性曲線濾波器的設計參數：通帶截止頻率fp=1300Hz;阻帶截止頻率 fs=1500Hz;（2）繪制出高通濾波器的幅頻和相頻特性曲線如下： 高通濾波器的幅頻和相頻特性曲線濾波器的設計參數：通帶截止頻率fp=300Hz;阻帶截止頻率 fs=100Hz;（3）濾波前信號的頻譜圖（4）經過低通濾波器后的語音信號頻譜圖 從圖像可以看出，原始聲音經過低通濾波器后，語音信號的幅度基本不變，濾波器把信號高頻部分濾掉了，只剩下低頻部分。用sound函數回放語音信號，發(fā)現原始聲音信號比濾波后聲音信號清晰，濾波后聲音顯得很低沉。（5）

11、經過高通濾波后的語音信號頻譜圖 從圖像可以看出，原始聲音經過高通濾波器后，語音信號的頻率基本不變，但一部分噪聲已經被濾掉，用sound函數回放語音信號，發(fā)現原始聲音信號與高通濾波后信號差別很小，但是聲音變得清晰，說明一部分噪聲已被濾掉。（6）經過帶通濾波器后的語音信號頻譜圖 從圖像可以看出，原始聲音經過帶通濾波器后，高頻部分被濾掉，部分噪聲被濾掉，用sound函數回放語音信號，原始聲音信號比濾波后聲音信號清晰，濾波后聲音比原始聲音低沉，但比經過低通濾波器后的清亮。（3） 求出特征頻段語音信號隨時間變化的曲線（每隔一定時間求一次功率譜，連接成曲線，即短時 FFT）。 我選取了語音信號出現較多的峰

12、值對應的頻率f=300Hz，500Hz，750Hz，900Hz，1000Hz，1200Hz，每隔0.1秒做功率譜，求出奇在017秒內的功率譜的變化情況，繪制出如下曲線：代碼：clcclear allFile=('I:Matlab大作業(yè)work.wav')Y1,Fs,nbits=wavread(File);Y1 = Y1(:,1)'%取雙聲道中的一個聲道l=length(Y1);Y=Y1(1:l); I = 300,500,750,900,1000,1200; %功率譜峰值位置p1 = zeros;p2 = zeros;p3 = zeros;p4 = zeros;p5 =

13、 zeros;p6 = zeros; for i = 1:170; S = Y(Fs/10*(i-1)+1 : Fs/10*i); z = fft(S,Fs); z1 = z(I(1); z2 = z(I(2); z3 = z(I(3); z4 = z(I(4); z5 = z(I(5); z6 = z(I(6); p1(i) = z1 .*conj(z1)/Fs; %功率譜 p2(i) = z2 .*conj(z2)/Fs; %功率譜 p3(i) = z3 .*conj(z3)/Fs; %功率譜 p4(i) = z4 .*conj(z4)/Fs; %功率譜 p5(i) = z5 .*conj

14、(z5)/Fs; %功率譜 p6(i) = z6 .*conj(z6)/Fs; %功率譜 endw = 0:17/(170-1):17; figure; plot(w,p1); xlabel('秒(s)'); title('300Hz功率譜變化圖像'); figure; plot(w,p2); xlabel('秒(s)'); title('500Hz功率譜變化圖像'); figure; plot(w,p3); xlabel('秒(s)'); title('750Hz功率譜變化圖像'); figure; plot(w,p4); xlabel('秒(s)'); title('900Hz功率譜變化圖像'); figure; plot(w,p5); xlabel('秒(s)'); title('1000Hz功率譜變化圖像'); figure; plot(w,p6); xlabel('秒(s)'); title('1