MATLAB语音FIR滤波器处理

%----------------------------------------------------------
%           MATLAB FIR语音滤波          
%
%         Copyright (c) 2011,Edward.xu        
%                             zhbit              
%----------------------------------------------------------
clear all;clc;
%语音采样
fs=8000;      %取样频率
duration=1;     %录音时间
fprintf('Press any key to start %g seconds of recording...\n',duration);
pause;
fprintf('Recording...\n');
data=wavrecord(duration*fs,fs);     %duration*fs 是总的采样点数
fprintf('Finished recording.\n');
fprintf('Press any key to play the recording...\n');
pause;
wavplay(data,fs);
%FFT变换
fs=8000;N=1024;  %采样频率和数据点数
n=0:N-1;t=n/fs;  %时间序列
x=data; %信号
y=fft(x,N);  %对信号进行快速Fourier变换
mag=abs(y);   %求得Fourier变换后的振幅
f=n*fs/N;  %频率序列
subplot(331);
plot(f,mag);  %绘出随频率变化的振幅
axis([0,4000,0,5]);
grid;
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');
title('原始信号FFT变换频谱');
%FIR低通滤波器
fp=1000;                 %通带截止频率
fs=1200;                 %阻带截止频率
At=100;                 %带外衰减系数(dB)
wave=1;                 %带内纹波系数(dB)
mval=[1 0];               %边界处的幅值
fedge=[fp fs];
derta1=1-10^(-wave/20);         %δ1的计算
derta2=10^(-At/20);           %δ2的计算
dev=[derta1 derta2];           %基本参数δ1与δ2
fs=8000;                 %FIR采样频率
%参数计算
[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs);
%雷米兹交替算法设计最优FIR滤波器
b=remez(N,fpts,mag,wt);
%FIR滤波器单位脉冲响应
[h,w]=freqz(b,1,1024);
subplot(332);
plot(w*(fs/2)/pi,20*log10(abs(h)));
grid;
xlabel('频率/Hz')
ylabel('幅度/dB')
title('FIR低通滤波器单位脉冲响应');
%FIR低通滤波器滤波
sf1=filter(b,1,data);
subplot(333);
fs=8000;N=1024;  %采样频率和数据点数
n=0:N-1;t=n/fs;  %时间序列
x=sf1; %信号
y=fft(x,N);  %对信号进行快速Fourier变换
mag=abs(y);   %求得Fourier变换后的振幅
f=n*fs/N;  %频率序列
plot(f,mag);  %绘出随频率变化的振幅
axis([0,4000,0,5]);
grid;
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');
title('FIR低通滤波后FFT变换频谱');
fprintf('Press any key to play the FIR LP filter recording...\n');
pause;
wavplay(sf1,fs);
%FIR高通滤波器单位脉冲响应
fp=2800;fs=3000;Rp=1;As=100;   % 滤波器指标
fb=[fp,fs];m=[0,1];        % 计算remezord函数所需参数f,m,dev
derta1=1-10^(-Rp/20);         %δ1的计算
derta2=10^(-As/20);           %δ2的计算
dev=[derta2 derta1];
[n,fo,mo,W]=remezord(fb,m,dev,8000);  % 确定remez函数所需参数
hn=remez(n,fo,mo,W);   % 调用remez函数进行设计,用于滤除噪声nt中的低频成分
%yt=filter(hn,1,1024);   %滤除随机噪声中低频成分，生成高通噪声yt
[h,w]=freqz(hn,1,1024);
subplot(334);
plot(w*(4000)/pi,20*log10(abs(h)));
grid;
xlabel('频率/Hz')
ylabel('幅度/dB')
title('FIR高通滤波器单位脉冲响应');
%FIR高通滤波器滤波
sf2=filter(hn,1,data);
subplot(335);
fs=8000;N=1024;  %采样频率和数据点数
n=0:N-1;t=n/fs;  %时间序列
x=sf2; %信号
y=fft(x,N);  %对信号进行快速Fourier变换
mag=abs(y);   %求得Fourier变换后的振幅
f=n*fs/N;  %频率序列
plot(f,mag);  %绘出随频率变化的振幅
axis([0,4000,0,5]);
grid;
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');
title('FIR高通滤波后FFT变换频谱');
fprintf('Press any key to play the FIR HP filter recording...\n');
pause;
wavplay(sf2,fs);
%FIR带通滤波器
fp1=1200;                 %通带截止频率
fc1=1000;                 %阻带截止频率
fp2=3000;
fc2=3200;
At=100;                 %带外衰减系数(dB)
wave=1;                %带内纹波系数(dB)
mval=[0 1 0];               %边界处的幅值
fedge=[fc1 fp1 fp2 fc2];
derta1=1-10^(-wave/20);         %δ1的计算
derta2=10^(-At/20);           %δ2的计算
dev=[derta2 derta1 derta2];       %基本参数δ1与δ2
fs=8000;                 %FIR采样频率
%参数计算
[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs);
%雷米兹交替算法设计最优FIR滤波器
b=remez(N,fpts,mag,wt);
%FIR滤波器单位脉冲响应
[h,w]=freqz(b,1,1024);
subplot(336);
plot(w*(fs/2)/pi,20*log10(abs(h)));
grid;
xlabel('频率/Hz')
ylabel('幅度/dB')
title('FIR带通滤波器单位脉冲响应');
%FIR带通滤波器滤波
sf3=filter(b,1,data);
subplot(337);
fs=8000;N=1024;  %采样频率和数据点数
n=0:N-1;t=n/fs;  %时间序列
x=sf3; %信号
y=fft(x,N);  %对信号进行快速Fourier变换
mag=abs(y);   %求得Fourier变换后的振幅
f=n*fs/N;  %频率序列
plot(f,mag);  %绘出随频率变化的振幅
axis([0,4000,0,5]);
grid;
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');
title('FIR带通滤波后FFT变换频谱');
fprintf('Press any key to play the FIR BP filter recording...\n');
pause;
wavplay(sf3,fs);