【正文】 頻譜： y_fil=filter(b,a,y)。 % 用設計好的濾波器對 y 進行濾波 Y_fil=abs(fft(y_fil))。Y_fil=Y_fil(1:N/2)。 % 計算頻譜取前一半 subplot(3,2,1)。plot(t,x)。xlabel(39。時間 (單位 :s)39。)。ylabel(39。幅度 39。)。title(39。原始吹管音樂 39。)。grid subplot(3,2,2)。plot(f,X)。xlabel(39。頻率 (單位 :Hz)39。)。ylabel(39。幅度譜 39。)。title(39。原始吹管音樂幅度譜圖 39。)。grid axis([0,4*10^3,0,2021]) subplot(3,2,3)。plot(t,y)。xlabel(39。時間 (單位 :s)39。)。ylabel(39。幅度 39。)。title(39。加噪后的吹管音樂39。)。grid subplot(3,2,4)。plot(f,Y)。xlabel(39。頻率 (單位 :Hz)39。)。ylabel(39。幅度譜 39。)。title(39。加噪后的吹管音鄒磊 《吹管音樂 濾波去噪 — 使用 雙線性變換設計切比雪夫 II型濾波器》 第 21頁 共 23頁 樂幅度譜圖 39。)。grid axis([0,4*10^3,0,2021]) subplot(3,2,5)。plot(t,y_fil)。xlabel(39。時間 (單位 :s)39。)。ylabel(39。幅度 39。)。title(39。濾波后的吹管音樂 39。)。grid axis([0,3,1,1]) subplot(3,2,6)。plot(f,Y_fil)。xlabel(39。頻率 (單位 :Hz)39。)。ylabel(39。幅度譜 39。)。title(39。濾波后的吹管音樂幅度譜圖 39。)。grid axis([0,4*10^3,0,2021]) figure subplot(211)。plot(f,20*log10(abs(Y)/max(abs(X))))。axis tight。grid on。 % 繪制加噪 吹管音樂 頻譜圖 subplot(212)。plot(f,20*log10(abs(Y_fil)/max(abs(X))))。axis tight。grid on。 % 繪制去噪 吹管音樂 頻譜圖 sound(y_fil,fs,bits)。 附錄二： afd_chb2 函數 function [b,a] = afd_chb2(Wp,Ws,Rp,As)。 % Analog Lowpass Filter Design: Chebyshev2 % % [b,a] = afd_chb2(Wp,Ws,Rp,As)。 % b = Numerator coefficients of Ha(s) % a = Denominator coefficients of Ha(s) % Wp = Passband edge frequency in rad/sec。 Wp 0 % Ws = Stopband edge frequency in rad/sec。 Ws Wp 0 % Rp = Passband ripple in +dB。 (Rp 0) % As = Stopband attenuation in +dB。 (As 0) % if Wp = 0 error(39。Passband edge must be larger than 039。) 鄒磊 《吹管音樂 濾波去噪 — 使用 雙線性變換設計切比雪夫 II型濾波器》 第 22頁 共 23頁 end if Ws = Wp error(39。Stopband edge must be larger than Passband edge39。) end if (Rp = 0) | (As 0) error(39。PB ripple and/or SB attenuation ust be larger than 039。) end ep = sqrt(10^(Rp/10)1)。 A = 10^(As/20)。 OmegaC = Wp。 OmegaR = Ws/Wp。 g = sqrt(A*A1)/ep。 N = ceil(log10(g+sqrt(g*g1))/log10(OmegaR+sqrt(OmegaR*OmegaR1)))。 fprintf(39。\n*** Chebyshev2 Filter Order = % \n39。,N) [b,a]=u_chb2ap(N,As,Ws)。 附錄三： 函數 function [db,mag,pha,grd,w] = freqz_m(b,a)。 % freqz 子程序的改進版本 % % [db,mag,pha,grd,w] = freqz_m(b,a)。 % db = [0 到 pi 弧度 ]區(qū)間內的相對振幅 (db) % mag = [0 到 pi 弧度 ]區(qū)間內的絕對振幅 % pha = [0 到 pi 弧度 ]區(qū)間內的相位響應 % grd = [0 到 pi 弧度 ]區(qū)間內的群遲延 % w = [0 到 pi 弧度 ]區(qū)間內的 501 個頻率樣本向量 % b = Ha(z)的分子多項式系數 (對 FIR b=h) % a = Ha(z)的分母多項式系數 (對 FIR: a=[1]) 鄒磊 《吹管音樂 濾波去噪 — 使用 雙線性變換設計切比雪夫 II型濾波器》 第 23頁 共 23頁 % [H,w] = freqz(b,a,1000,39。whole39。)。 H = (H(1:1:501))39。 w = (w(1:1:501))39。 mag = abs(H)。 db = 20*log10((mag+eps)/max(mag))。 pha = angle(H)。 % pha = unwrap(angle(H))。 grd = grpdelay(b,a,w)。 % grd = diff(pha)。 % grd = [grd(1) grd]。 % grd = [0 grd(1:1:500)。 grd。 grd(2:1:501) 0]。 % grd = median(grd)*500/pi。