[x1,fs]=audioread('G:\G盘恢复数据备份,2018.9.17\2222222222222222222222222222222222\Root\1.刘金豆\留下\个人(这里面有些东西时间有些久远了,部分是半成品,部分文件里面很乱,我也不太能分清了,大家取有用的就好)\课程设计\大三上学期课设1\1yuanshi.m4a');
f=fs*(0:32767)/8182;
X1=fft(x1,65536);w=(0:32767)/32768*20000;
%sound(x1,16000);
figure(1)
subplot(2,1,1),plot(x1,'b');title('FIR 低通滤波器滤波前的时域波形x1');
subplot(2,1,2),plot(w,abs([X1(1:32768)]));xlabel('频率/Hz');ylabel('幅值');
title('FIR 低通滤波器滤波前的频谱')
%窗函数法设计FIR数字滤波器
Ts=1/fs;
wp=2*pi*2000/fs;
ws=2*pi*2200/fs;
Rp=3;
Rs=55;
wdelta=ws-wp;
N=ceil(8*pi/wdelta);wn=(ws+wp)/2;
[b,a]=fir1(N,wn/pi,hamming(N+1));
figure(2)
freqz(b,a,5120);
title('FIR 低通滤波器');
x2=filter(b,a,x1);
X2=fft(x2,65536);w=(0:32767)/32768*20000;
sound(x2,16000);
figure(3)
subplot(2,1,1),plot(x2,'b');title('FIR 低通滤波器滤波后的时域波形x2');
subplot(2,1,2),plot(w,abs([X2(1:32768)]));xlabel('频率/Hz');ylabel('幅值');
title('FIR 低通滤波器滤波后的频谱')
xO=x1-x2;
XO=fft(xO,65536);w=(0:32767)/32768*20000;
figure(4)
subplot(2,1,1),plot(xO,'b');title('波形差x1-x2');
subplot(2,1,2),plot(w,abs([XO(1:32768)]));xlabel('频率/Hz');ylabel('幅值');
title('波形差的频谱')