[x,fs,nbit]=wavread('Test.wav');%讀取Test音頻文件
sound(x,fs)
N=length(x);
n=1:N;
c=0.1*sin(2*pi*1000*n./fs);
x1=x'+c;%將正弦波加載于原音頻文件
sound(x1,fs);
subplot(2,2,1);stem(n,x);
xlabel('采樣點');
ylabel('模擬輸入 (V)');
title('原始聲音信號'),xlabel('n')
X=fft(x);magX=abs(X(1:(N+1)/2+1));
k=0:1:(N+1)/2;w=2*pi/N*k;
subplot(2,2,2);stem(w/pi,magX);
title('原始聲音信號幅頻圖');xlabel('以pi為單位的頻率')
subplot(2,2,3);stem(n,x1);
xlabel('采樣點');
ylabel('模擬輸入 (V)');
title('加入噪聲的聲音信號');xlabel('n')
X1=fft(x1);magX1=abs(X1(1:(N+1)/2+1));
subplot(2,2,4);stem(w/pi,magX1);
title('加入噪聲信號幅頻圖');xlabel('以pi為單位的頻率')
wp=0.02*pi;%通帶截止頻率
ws=0.04*pi;%阻帶截止頻率
Rp=1;%通帶紋波
As=20;%阻帶紋波
T=1;
OmegaP=(2/T)*tan(wp/2);%原型通帶頻率
OmegaS=(2/T)*tan(ws/2);%原型阻帶頻率
ep=sqrt(10^(Rp/10)-1);%通帶紋波參數(shù)
Ripple=sqrt(1/(1+ep*ep));%通帶紋波
Attn=1/(10^(As/20));%阻帶衰減
[cs,ds]= p_butt(OmegaP,OmegaS,Rp,As);
[b,a]= bilinear(cs,ds,T);%雙線性變換
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))'; w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);
y=filter(b,a,x1);%b,a為巴特沃茲低通濾波器分子、分母系數(shù)
Y1=fft(y);magY1=abs(Y1(1:(N+1)/2+1));k=0:1:(N+1)/2;w=2*pi/N*k;
stem(w/pi,magY1);title('帶噪聲信號的濾波后的幅頻圖')
sound(y,fs);