2014用巴特莱特窗函数法设计数字FIR低通滤波器

2.FIR数字滤波器设计的基本步骤如下：

（1）确定技术指标

在设计一个滤波器之前,必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中,数字滤波器常被用来实现选频操作。因此,指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以2种方式给出。第一种是绝对指标。他提供对幅度响应函数的要求,一般应用于FIR滤波器的设计。第二种指标是相对指标。他以分贝值的形式给出要求。本文中滤波器的设计就以线性相位FIR滤波器的设计为例。

（2）逼近

确定了技术指标后,就可以建立一个目标的数字滤波器模型（通常采用理想的数字滤波器模型）。之后,利用数字滤波器的设计方法（窗函数法、频率采样法等）,设计出一个实际滤波器模型来逼近给定的目标。

（3）性能分析和计算机仿真

上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性,以验证设计结果是否满足指标要求;或者利用计算机仿真实现设计的滤波器,再分析滤波结果来判断。

三 程序流程图

开始 ↓

读入窗口长度 ↓ 计算hd(n)( ↓

调用窗函数子程序找w(n) ↓

计算h(n)=hd(n)w(n) ↓

调用子程序计算H（k）=DFT[h(n)] ↓

调用绘图子程序绘制H（k）幅度相位曲线 ↓

结束

四 程序源代码

%实验设计程序如下：

clc;

fp=500;fs=400;Fs=2000; Rp = 1; Rs = 40;

wp=2*pi*fp/Fs;ws=2*pi*fs/Fs; i=1-10^(-Rp/20); j=10^(-Rs/20);

%N=kaiord(0.1087, 0.0100,500,400,2000)

n=-20*log(sqrt(i*j))-13; N=n/14.6/(wp-ws)*2*pi

floor(N); N=76;

wn=(wp+ws)/2/pi;

b=fir1(N,wn,bartlett(N+1)); figure(1)

[h,w]=freqz(b,1,512,2000); g = 20*log10(abs(h));

plot(w,g);grid on; axis([0 1000 -50 3]);

xlabel('频率，Hz'); ylabel('增益，title('巴特莱特 LPF'); figure(2) zplane(b,1);

xlabel('b'); ylabel('a'); title('传输零极点');

f1=300;f2=600; n=0:600; t=n/10000;

x1=sin(2*pi*f1*t); x2=sin(2*pi*f2*t); x=x1+x2; figure(3)

subplot(2,1,1); plot(x1);grid on;

axis([0,50*pi,-3,3]);

xlabel('t');ylabel('x1'); title('x1的波形'); subplot(2,1,2); plot(x2);grid on;

axis([0,50*pi,-3,3]);

xlabel('t');ylabel('x2'); title('x2的波形'); figure(4)

subplot(2,1,1); plot(x);grid on;

axis([0,50*pi,-3,3]); xlabel('t');ylabel('x'); title('输入x的波形'); y=filter(b,1,x); subplot(2,1,2) plot(y);grid on;

) dB' axis([0,50*pi,-5,5]); xlabel('t');ylabel('y'); title('滤波器输出y的波形'); %频谱图

fs=2000;N=1024; n=0:N-1;t=n/fs; Y1=fft(x,N); Y2=fft(y,N);

mag1=abs(Y1);mag2=abs(Y2); f=n*fs/N; figure(5)

subplot(2,1,1);

plot(f(1:N/2),mag1(1:N/2)); title('输入信号的频谱图');

xlabel('频率/HZ');ylabel('振幅'); grid on;

subplot(2,1,2);

plot(f(1:N/2),mag2(1:N/2)); title('输出信号的频谱图');

xlabel('频率/HZ');ylabel('振幅'); grid on;

function N = kaiord(Fp, Fs, Rp, Rs, FT) % Computation of the length of a linear-phase % FIR multi-band filter using Kaiser's formula %

% Rp is the passband ripple % Rs is the stopband ripple % Fp is the passband edge in Hz % Fs is the stopband edge in Hz % FT is the sampling frequency in Hz. % If none specified default value is 2 % N is the estimated FIR filter order if nargin == 4, FT = 2; end

if length(Fp) > 1,

TBW = min(abs(Fp(1) - Fs(1)), abs(Fp(2) - Fs(2))); else

TBW = abs(Fp - Fs); end

num = -20*log10(sqrt(Rp*Rs))-13; den = 14.6*TBW/FT; N = ceil(num/den);