电子系统课程设计报告

FIR数字滤波器设计

专    业：   07级电子信息工程      

班    级：       电信一班      

学    号：     XXXXXXX     

姓    名：        XXX       

指导老师：    XXXXXXXXX    

同组成员：        XXX        

时间：20##年7月1日

一 、设计目的

1、熟悉FIR滤波器设计的基本方法。

2、掌握用窗函数设计FIR滤波器的基本原理和方法，熟悉MATLAB语言。

3、熟悉线性相位FIR滤波器的幅频特性和相位特性。

4、了解各种不同窗函数对滤波器性能的影响。

5、熟悉MCS51系列单片机基本原理。

6、掌握MCS51系列单片机基本指令和编程。

7、利用MATLAB仿真的数据，用单片机实现FIR滤波器的设计。

二、原理和方法

（一）FIR滤波器的设计

滤波器是设计来进行频率选择或频率分辨任务的线性时不变系统的通用名称。离散时间LTI(线性时不变)系统称为数字滤波器。按照LTI系统的单位脉冲响应长度是否有限可以分为有限长度脉冲响应(FIR)滤波器和无限长度脉冲响应(IIR)滤波器。

FIR滤波器可以利用快速傅立叶变换进行快速卷积来实现，另外，FIR滤波器具有严格的相位特性，这对于语音信号处理和数据传输是很重要的。目前FIR滤波器的设计方法主要有三种：窗函数法、频率采样法和切比雪夫等波纹逼近的最优化设计方法。最常用的是窗函数法和切比雪夫等波纹逼近的最优化设计法。本设计中用的窗函数法比较简单，可以应用现成的窗函数公式，在技术要求不高的时候是比较方便灵活的。它从时域出发，用一个窗函数截取理想的hd(n)得到h(n)，以有限长序列h(n)近似理想的hd(n)；如果从频域出发，用理想的Hd(ejw)在单位圆上等角度取样得到h(k)，根据h(k)得到h(z)将逼近理想的hd(z)，这就是频率取样法。

（二）窗函数法原理

众所周知，数字滤波器的设计是基于一定的性能指标。窗函数法设计滤波器也不例外，首先要对滤波器提出性能指标。一般是给定一个理想的频率响应Hd(ejw)，使所设计的FIR滤波器的频率响应H(ejw)去逼近所要求的理想的滤波器的响应Hd(ejw)。窗函数法设计的任务在于寻找一个可实现（有限长单位脉冲响应）的传递函数：

                     （4）

去逼近Hd(ejw)。一般情况下，一个理想的频率响应Hd(ejw)的傅立叶反变换：

                   （5）

所得到的理想单位脉冲响应hd(n)是一个无限长序列。对hd(n)经过适当的加权、截断处理才能得到一个所需要的有限长脉冲响应序列。对应不同的加权、截断，就有不同的窗函数。而所要寻找的滤波器脉冲响应就等于理想脉冲响应和窗函数的乘积，即：

                         （6）

由此可见，窗函数的性质就决定了滤波器的品质。由于窗函数有一个等于M的有限长度，所以它的响应有一个主瓣，其宽度与1/M成正比，还有高度比较小的旁瓣。主瓣在Hd(ejw)中产生一个过渡带，主瓣的宽度决定过渡带的宽度，即与1/M成正比，主瓣越宽，则过渡带越宽，旁瓣产生波动。

以下是几种常用的窗函数：

1、矩形窗

                   （7）

    2、汉宁(Hanning)窗

               （8）

   3、哈明(Hamming)窗

               （9）

   4、布莱克曼(Blackman)窗

        （10）

   5、凯泽(Kaiser)窗

           （11）

其中I0[.]是修正过的零阶贝塞尔(Bessel)函数，β为依赖于M的参数，选择M可以产生各种过渡带宽和接近最优的阻带衰减。对于相同的M，凯泽窗可以提供不同的过渡带宽，这是其它窗所做不到的。

窗函数法设计线性相位FIR滤波器可以按如下步骤进行：

（1）确定数字滤波器的性能要求。确定各临界频率{ωk}和滤波器单位脉冲长度M。

（2）根据性能要求和M值，合理地选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应Hd(ejw)的幅频特性和相位特性。

（3）用傅立叶反变换公式（5）求得理想单位脉冲响应hd(n)。

（4）选择适当的窗函数w(n)，根据式（6），求得所设计的FIR滤波器单位脉冲响应。

（5）用傅立叶变换求得其频率响应H(ejw)，分析它的幅频特性，若不满足要求，可适当改变窗函数形式或长度M，重复上述过程，直到得到满意的结果。

三、MATLAB实现的内容及步骤

（一）编制实验用主程序及相应子程序

1、在实验编程之前，认真复习有关FIR滤波器设计的有关知识，尤其是窗函数法的有关内容，阅读本实验原理与方法，熟悉窗函数及四种线形相位FIR滤波器的特性，掌握窗函数设计滤波器的具体步骤。

2、编制窗函数设计FIR滤波器的主程序及相应子程序。

（1）根据Hd(ejw)，利用傅立叶变换，编制产生理想单位脉冲响应hd(n)的简单子程序。

（2）根据公式（7）到公式（10）编制相应的窗函数序列的子程序。

1＞矩形窗程序：

  function w=jvxing(x);

M=x;

for n=1:M;

      w(n)=1;

end

2＞汉宁窗程序：

function w=hn(x)

M=x

for n=1:M

      w(n)=0.5*[1-cos(2*pi*n/(M+1))]

end

3＞汉明窗程序：

function w=hm(x);

M=x

for n=1:M

      w(n)=0.54-0.46*cos(2*pi*(n-1)/(M-1))

end

4＞布莱克曼床程序：

function w=blkm(x);

M=x

for n=1:M

      w(n)=0.42-0.5*cos(2*pi*(n-1)/(M-1))+0.08*cos((4*pi*(n-1)/(M-1)))

end

5＞凯泽窗程序：

function w=kz(x);

M=x;

beta=4（6，8）;

for n=1:M;

      w(n)=besseli(0,(beta*sqrt(1-((2*n/(M-1))-1)^2)))/besseli(0,beta);

end

窗函数图形：

（3）在上述子程序的基础上设计主程序完成线性相位  FIR滤波器的窗函数法设计。

（二）上机实验内容

1、用Hanning窗设计一个线形相位带通滤波器，其长度M=30，技术指标要求如下：

低端阻带边缘：ω1s=0.2π ;     低端通带边缘：ω1p=0.35π

高端通带边缘：ω1p=0. 65π  ;     高端阻带边缘：ω1s=0.8π

绘制幅频和相位特性曲线、理想脉冲响应和实际脉冲响应离散数据图。如果M=45，重复这一设计，观察幅频和相位特性的变化，注意长度M变化对结果的影响。

主程序：

Wp1=0.35*pi;Wp2=0.65*pi;Ws1=0.2*pi;Ws2=0.8*pi;

M=30;

N=M-1;

M1=45;

N1=M1-1;

w=hn(M);                                   %汉宁窗w1=hn(M1);

Wc1=(Wp1+Ws1)/2;Wc2=(Wp2+Ws2)/2;

k=0:N;

k1=0:N1;

hd=(Wc2/pi)*sinc(Wc2*(k-0.5*N)/pi)-(Wc1/pi)*sinc(Wc1*(k-0.5*N)/pi);

hd1=(Wc2/pi)*sinc(Wc2*(k1-0.5*N1)/pi)-(Wc1/pi)*sinc(Wc1*(k1-0.5*N1)/pi);

h=hd.*w;  h1=hd1.*w1;

omega=linspace(0,pi,512);

mag=freqz(h,1,omega);magdb=20*log(abs(mag));

mag1=freqz(h1,1,omega);magdb1=20*log(abs(mag1));

subplot(321)

plot(omega/pi,magdb,'-b',omega/pi,magdb1,'-r');

title('幅频特性（加汉宁窗）')

xlabel('normalizek frequency (rad/s)*π');

ylabel('Gain,dB');

legend('M=30','M=45');

grid;

phase=unwrap(angle(mag));

phase1=unwrap(angle(mag1));

subplot(322);

plot(omega,phase,'-b',omega,phase1,'-r');

title('相位特性');

xlabel('Normalized frequency rad/s');

ylabel('phase rad');

legend('M=30','M=45');grid;

subplot(323)

stem(k,h);

title('实际脉冲响应离散数据图M=30');

grid;

subplot(324)

stem(k,hd);

title('理想脉冲响应离散数据图M=30');

grid;

subplot(325)

stem(k1,h1);

title('实际脉冲响应离散数据图M=45');

grid;

subplot(326)

stem(k1,hd1);

title('理想脉冲响应离散数据图M=45');

grid;

图形：

分析：随着长度M的增加滤波器的过渡带减小，阻带特性变好，理想脉冲和实际脉冲幅度增大。

2、改用矩形窗和Blackman窗，设计步骤（1）中的带通滤波器，观察并记录窗函数对滤波器幅频和相位特性的影响，比较这三种窗函数的特点。

程序：

Wp1=0.35*pi;Wp2=0.65*pi;Ws1=0.2*pi;Ws2=0.8*pi;

N=45;

M=N-1;

w1=jvxing(N);

w2=hn(N);

w3=blkm(N);

Wc1=(Wp1+Ws1)/2;Wc2=(Wp2+Ws2)/2;

k=0:M;

hd=(Wc2/pi)*sinc(Wc2*(k-0.5*M)/pi)-(Wc1/pi)*sinc(Wc1*(k-0.5*M)/pi);

h1=hd.*w1; 

h2=hd.*w2;

h3=hd.*w3; 

omega=linspace(0,pi,512);

mag1=freqz(h1,1,omega);magdb1=20*log(abs(mag1));

mag2=freqz(h2,1,omega);magdb2=20*log(abs(mag2));

mag3=freqz(h3,1,omega);magdb3=20*log(abs(mag3));

subplot(211)

plot(omega/pi,magdb1,'b-',omega/pi,magdb2,'r-',omega/pi,magdb3,'g-');

legend('矩形窗','汉纳窗','布莱克曼窗')

title('幅频特性')

xlabel('normalizek frequency');

ylabel('Gain,dB');

grid;

phase1=unwrap(angle(mag1));

phase2=unwrap(angle(mag2));

phase3=unwrap(angle(mag3));

subplot(212)

plot(omega,phase1,'b-',omega,phase2,'r-',omega,phase3,'g-');

legend('矩形窗','汉纳窗','布莱克曼窗')

title('相位特性');

ylabel('phase');

grid;

图形：

分析：在长度相同（M=45）情况下，加矩形窗的滤波器过渡带要小于加布莱克曼窗，而加布莱克曼窗的阻带特性比加矩形窗好（衰减大），加汉纳窗则介于二者之间。

相位特性：加矩形窗、汉纳窗、布莱克曼窗滤波器在通带相位特性良好（基本是直线）。

3、用Kaiser窗设计一个专用的线性相位滤波器。M=40，理想幅频特性如下：

当β值分别为4,6,8时，设计相应的滤波器，比较它们的幅频和相位特性，观察并分析β值不同的时候对结果有什么影响。

主程序：

M=40;

Wc1=0.2*pi;Wc2=0.4*pi;Wc3=0.6*pi;Wc4=0.8*pi;

w=kz(M);                                             %凯泽窗

w1=kz1(M);

w2=kz2(M);

k=1:M;

hd=(Wc2/pi)*sinc(Wc2*(k-0.5*M)/pi)-(Wc1/pi)*sinc(Wc1*(k-0.5*M)/pi)+(Wc4/pi)*sinc(Wc4*(k-0.5*M)/pi)-(Wc3/pi)*sinc(Wc3*(k-0.5*M)/pi);

h=hd.*w;

h1=hd.*w1;

h2=hd.*w2;

omega=linspace(0,pi,512);

mag=freqz(h,1,omega);                                  %求h的频率响应

mag1=freqz(h1,1,omega);

mag2=freqz(h2,1,omega);

magdb=20*log10(abs(mag));

magdb1=20*log10(abs(mag1));

magdb2=20*log10(abs(mag2));

subplot(121)

plot(omega/pi,magdb,'-b',omega/pi,magdb1,'-r',omega/pi,magdb2,'-g');

title('频率特性')

xlabel('Normalized frequency (rad/s)*π');

ylabel('Gain dB');

legend('β=4','β=6','β=8');

grid;

phase=unwrap(angle(mag));                              %求mag的相位

phase1=unwrap(angle(mag1));

phase2=unwrap(angle(mag2));

subplot(122);

plot(omega,phase,'-b',omega,phase1,'-r',omega,phase2,'-g');

title('相位特性');

xlabel('Normalized frequency rad/s');

ylabel('phase rad');

legend('β=4','β=6','β=8');

grid;

图形：

分析：随着β的增大，滤波器两边的阻带衰减增大（阻带特性变好），过渡带增大，导致中间阻带特性变差。通带传递函数的相位良则不受β变化的影响。

四、DVCC系列仿真系统实现滤波器

程序流程图：

主程序：

ORG 0200H

LJMP START

TABLEHX:DB 05h,04h,03h,02h,01h,00h,00h,00h

TABLEX: DB 01h,03h,01h,02h,00h,00h,00h,00h

TABLE:DB 0C0H,0F9H,0A4H,0B0H,99H,92H,82H,0F8H,80H,90H

ORG 0500H

START:  MOV R0,#30H

        MOV R2,#0H

        MOV R3,#8H

s:     MOV A,R2

    MOV DPTR,#TABLEX

    MOVC A,@A+DPTR

    MOV @R0,A

        INC R0

    INC R2

    DJNZ R3,s

      MOV R1,#40H  ;存Y(N)

      MOV R3,#8H   ;乘的次数

          MOV R4,#0H       ;N

s1:       MOV R6,#0H      ;Y(N)=x[k]h[n-k]

          MOV DPTR,#TABLEHX

      MOV A,R4

      ADD A,#01H

      MOV R5,A

          MOV R2,#0H     ;取H(↘)即K

s2:       MOV A,R2

          MOVC A,@A+DPTR    ;从H取出

          MOV R7,A      ;暂存H(K)

          MOV A,R4       ;赋N--＞A

          CLR C

          SUBB A,R2      ;N-K

          ADD A,#30H

          MOV R0,A

          MOV A,@R0    ;取X(N-K)

          MOV B,R7      ;暂存H[K]

          MUL AB

          ADD A,R6       ;加上[K]*X[N-K]

          MOV R6,A        ;暂存H[K]*X[N-K]

          INC R2      ;K

          DJNZ R5,s2   ;R5:要加的乘积项数

          MOV A,R6

          MOV @R1,A     ; 存Y[N]

      INC R1         ; Y[N+1]

      INC R4         ;N

          DJNZ R3,s1   ;R3乘的次数

        ;;;;;;;;;;;;;;

       

dx1:    MOV R6,#0H      ;N循环置初

dx:     MOV A,#40H

        ADD A,R6

        MOV R0,A

        MOV A,@R0     ;取Y[]

        MOV B,#10

        DIV AB

        MOV 7DH,B       ;个位，可在范围内任取

        MOV B,#10

        DIV AB

        MOV 7EH,B       ;十

        MOV 7FH,A       ;百

        MOV R0,#21H      ; 数码管选通端口地址

        MOV R1,#22H      ;数据输出端口地址

        MOV DPTR,#TABLE     ;段位码

LED:    MOV A,7DH

        MOVC A,@A+DPTR

        MOVX @R1,A

        MOV A,#01H      ;  选数码管 左6

        MOVX @R0,A     ;显示个位

        ACALL DELAY

         MOV A,7EH

        MOVC A,@A+DPTR

        MOVX @R1,A

        MOV A,#02H       ;选数码管 左5

        MOVX @R0,A       ;显示十位

        ACALL DELAY

         MOV A,7FH

        MOVC A,@A+DPTR

        MOVX @R1,A

        MOV A,#04H        ;选数码管 左4

        MOVX @R0,A        ;显示百位

        ACALL DELAY

        ;;;;;;;;;;;;;

        MOV 6AH,R6         ;第幾位

        MOV A,6AH

        MOV B,#0AH          ;B中放10

        DIV AB

        MOV 68H,B

        MOV 69H,A

        MOV R7,#0FFH

       

        MOV A,68H

        MOVC A,@A+DPTR

        MOVX @R1,A

        MOV A,#10H           ;选数码管左2

        MOVX @R0,A           ;显示个数

        ACALL DELAY

        MOV A,69H

        MOVC A,@A+DPTR

        MOVX @R1,A

        MOV A,#20H           ;选数码管左1

        MOVX @R0,A           ;显示十位

        ACALL DELAY

        ;;;;;;;;;;;;;;;;

        ACALL KEY1

        CJNE R6,#8H,dx       ;R6乘的次数

        LJMP dx1             ;对R6清零 ,使7变为0

KEY1:   ACALL KS1

        JNZ LK1

    RET

LK1:     ACALL DELAY3

        INC R6

        RET

KS1:    MOV DPTR,#0FF21H             

    MOV A,#00H

        MOVX @DPTR,A

        MOV DPTR,#0FF23H             

        MOVX A,@DPTR

    CPL A

    ANL A,#0FH          ;查看是否有按键       

    RET

DELAY:

        MOV R5,#02H

D1:     MOV R4,#0FFH

        DJNZ R4,$

        DJNZ R5,D1

        RET

DELAY3:  MOV R7,#07FH

DD1:      ACALL DELAY

        DJNZ R7,DD1

        RET

END

运行结果：00 005，01 019，02 020，03 025，

          04 018，05 011，06 005，07 002

五、思考题

1、如果没有给定h(n)的长度M，而是给定了带通边缘截止频率ωc和阻带临界频率ωs，以及相应的衰减，你能根据这些条件用窗函数法设计线性相位FIR低通滤波器吗？

答：在给定带通边缘截止频率ωc和阻带临界频率ωs情况下，可以通过查下表，依据不同的窗函数，近似求得长度M，再由相应的衰减就可得出FIR低通滤波器。

   表（1） 常用窗函数特性

2、如果在运用（5）式产生单位脉冲响应h(n)时Hd(ejw)比较复杂，怎样才能产生实际脉冲响应h(n)？

答：可以用频率取样法。

对Hd(ejΩ)进行频率抽样得到Hd[m]的值，再对Hd[m]做M+1点IDFT,就可以得到实际脉冲响应h(n)

3、频率取样法和窗函数法各有什么特点？（简单说明）

答：频率采样法对目标系统在某些频率点ωk进行取样，要求所设计的FIR滤波器的频率响应在这些指定的频率点ωk上的频率特性与目标系统在这些点上的特性完全相同，以此为依据设计出合适的FIR滤波器。

窗函数法依据给出的Hd(ejw)函数，通过窗函数的截断，从而得到H(ejw)函数，根据这个函数，可以通过傅里叶反变换直接计算系统的单位响应h(n)并由此构建数值滤波器。窗函数设计法的基本思想，是找出一个频率特性为Hd(ejw)的、与H(ejw)尽可能接近的FIR函数。窗函数的性质就决定了滤波器的品质。由于窗函数有一个等于M的有限长度，所以它的响应有一个主瓣，其宽度与1/M成正比，还有高度比较小的旁瓣。主瓣在Hd(ejw)中产生一个过渡带，主瓣的宽度决定过渡带的宽度，即与1/M成正比，主瓣越宽，则过渡带越宽，旁瓣产生波动

4、频率特性图中的横轴单位是什么？

答：横轴单位是（rad/s）*π