数字信号处理之低通滤波器设计

 1 clc;
 2 clear;
 3 close all;
 4   
 5 [x1,fs,bits]=wavread('qq.wav'); %读取语音信号的数据，赋给变量x1
 6 L=length(x1);
 7 sound(x1,fs); %播放语音信号
 8 
 9 N=L;
10 k=[0:N-1];
11 y1=fft(x1,N);  %对信号做L个点FFT变换
12 figure(1);
13 subplot(2,1,1);
14 plot(k(1:20000),abs(y1(1:20000))); %做原始语音信号的FFT频谱图
15 title('原始语音信号FFT频谱');
16   
17 f=fs*(0:N-1)/N;
18 subplot(2,1,2);
19 plot(f(1:20000),abs(y1(1:20000)));
20 title('原始语音信号频谱');
21 xlabel('f Hz');ylabel('fuzhi');

 1 %给原始的语音信号加上一个高频余弦噪声，频率为5kHz。画出加噪后的语音信号时域和频谱图，与原始信号对比，可以很明显的看出区别。
 2 t=[0:1/fs:(L-1)/fs];        %将所加噪声信号的点数调整到与原始信号相同，构造采样时间点（模拟时间）
 3 Au=0.03;                    %噪声幅度
 4 d=[Au*cos(2*pi*8000*t)]';   %噪声为8kHz（2*pi*8000/(2*pi)）的余弦信号(模拟时间)
 5 x2=x1+d;
 6 
 7 figure(2);
 8 subplot(2,1,1);
 9 plot(x1);        %做原始语音信号的时域图形
10 title('原始语音信号');xlabel('time n');ylabel('fuzhi');
11 
12 sound(x2,fs);    %播放加噪声后的语音信号
13 subplot(2,1,2);
14 plot(x2);        %做原始语音信号的时域图形
15 title('加噪后的信号');xlabel('time n');ylabel('fuzhi');
16 
17 y2=fft(x2,N);     %对信号做L个点FFT变换
18 figure(3);
19 subplot(2,1,1);
20 plot(k(1:20000),abs(y2(1:20000))); %做原始语音信号的FFT频谱图
21 title('加躁语音信号FFT频谱');
22   
23 subplot(2,1,2);
24 plot(f(1:20000),abs(y2(1:20000)));
25 title('加躁语音信号频谱');
26 xlabel('f Hz');ylabel('fuzhi');

 1 wp=2*pi*4000;                       %通带边界角频率
 2 ws=2*pi*5000;                       %阻带边界角频率
 3 Rp=1;                               %通带最大衰减
 4 Rs=15;                              %阻带最小衰减
 5 [NN,Wn]=buttord(wp,ws,Rp,Rs,'s');   %选择滤波器的最小阶数
 6 [Z,P,K]=buttap(NN);                 %创建butterworth模拟滤波器
 7 [Bap,Aap]=zp2tf(Z,P,K);
 8 [b,a]=lp2lp(Bap,Aap,Wn);   
 9 [bz,az]=bilinear(b,a,fs);          %用双线性变换法实现模拟滤波器到数字滤波器的转换
10 
11 [H,W]=freqz(bz,az);                %绘制频率响应曲线
12 figure(4);
13 plot(W*fs/(2*pi),abs(H));
14 grid;
15 xlabel('频率／Hz');
16 ylabel('频率响应幅度');
17 title('Butterworth');

x3=filter(bz,az,x2);
figure(5);
subplot(2,1,1);
plot(t,x2);               %画出滤波前的时域图
title('滤波前的时域波形');
subplot(2,1,2);
plot(t,x3);               %画出滤波后的时域图
title('滤波后的时域波形');
sound(x3,fs);        %播放滤波后的信号

y2=fft(x2,N);  %对信号做L个点FFT变换
figure(6);  
subplot(2,1,1);
plot(f(1:20000),abs(y2(1:20000)));
title('加躁语音信号频谱');
xlabel('f Hz');ylabel('fuzhi');

y2=fft(x3,N);  %对信号做L个点FFT变换
subplot(2,1,2);
plot(f(1:20000),abs(y2(1:20000)));
title('去噪语音信号频谱');
xlabel('f Hz');ylabel('fuzhi');