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);