%巴特沃斯濾波器
clear all;
close all;
fp=1000;%通帶截止頻率
fs=2000;%阻帶截止頻率
Rp=2;
Rs=35;
Fs=8000;%設(shè)定采樣頻率
Wp=2*fp/Fs;
Ws=2*fs/Fs;
[N,Wn]=buttord(Wp,Ws,Rp,Rs);
[b,a]=butter(N,Wn);
[H,w]=freqz(b,a,Fs);
h=abs(H);
plot(w/pi,h);
axis([0 5 0 1.1]);
xlabel('w/pi');
ylabel('|H(jw)|');
title('巴特沃斯濾波器');
grid on
 
%低頻發(fā)生器產(chǎn)生原信號
t=[0:100]/Fs;
[x Fs Nbits]=wavread('myvoice1');
%A=1;
%f=1000;
%x=A*sin(2*pi*f*t);
figure(2);
plot(x(1:100)); 
title('原信號');
aver_x=mean(x)
power_x=mean(x.^2)
d_x=var(x)
xf=fft(x,Fs);
figure(3);
stem(abs(xf),'.')
title('原信號的頻譜');
rxx=xcorr(x,'coef');
figure(4)
plot(rxx)
title('原信號的自相關(guān)')
rxxf=fft(rxx,Fs);
figure(5);
stem(abs(rxxf),'.');
title('原信號功率譜密度');
 
 
%輸入信號,加上噪聲
xn=AWGN(x,10,0);  
figure(6);
plot(t,xn(1:length(t)));
title('加噪后信號');
 

%通過低通濾波器
y1=filter(b,a,xn);
figure(7);
plot(t,y1(1:length(t)));
title('濾波器輸出信號');
figure(8)
for i=1:21600
    y(i)=mean(y1(i:8:21600));
end
plot(y(1:100))
title('提取信號')
aver_y=mean(y)
power_y=mean(y.^2)
d_y=var(y)
yf=fft(y,Fs);
figure(9);
stem(abs(yf),'.');
title('輸出信號的頻譜');
ryy=xcorr(y,'coef');
ryyf=fft(ryy,Fs);
figure(10);
stem(abs(ryyf),'.');
title('輸出信號功率譜密度');
rxy=xcorr(x,y,'coef');
figure(11);
stem(rxy,'.');
title('互相關(guān)函數(shù)');