現代信號處理大型作業

一.試用奇階互補法設計兩帶濾波器組(高、低通互補)，進而實現四帶濾

波器組；并畫出其頻響。濾波器設計參數為：Fp=1.7KHz,F「=23KHz,Fs

=8KHz,加侖70dB0

(一)、分析

與通常的濾波器相比，互補濾波器具有優良的結構特性和結構特性，具

有較低的噪聲能量和系數敏感性，其定義如下：

一組濾波器H1(Z),42(Z),…….H.(Z)如果滿足下式:

Z|"K(〃W)|=1,0<w<2兀

則稱這組濾波器為幅度互補濾波器;如果滿足下式:

=1,0<W<2K

則稱這組濾波器為功率互補濾波器，同時互補濾波器還應該滿足：

XHk(z)=A(z)其中A(z)為全通函數，適當的選擇全通函數，可以

使兩帶函數具有所需要的低通和高通特性。

(二)、設計步驟

(1)對Fp、Fr進行預畸

V\Fp

=爾（

Q=應（

(2)計算。：=叵而7,判斷Q；是否等于1，即該互補濾波器是否為互補

鏡像濾波器

(3)計算相關系數

jll=[(10°Upmin-1)(10°Mrmin-1)",]2；

,Qp

k="一人

ii-VF

q=q°+2加+15/+150q；;

N=lg(42/i6)/lgq；

i；(N活數)

〃[，T(N為禺數)

2公￡(_1尸產7sin((2m+l>f/)

，"=oN

l+2X(-ir^2COS(2〃77r

m=lN

v,=”-吟)(1-C〃Q;

a=1,…M

i=i+c"

夕l&E…M

(4)互補鏡像濾波器的數字實現

2-a,

A=

2+a.

2-A

B

i=2+0

口4+z々

H.(Z)=1,…M

1J1+AZ-2

pI7-2

區“心口七毛八卜以

j1TUj乙

//JZ)=1[/7.(Z)+/72(Z)];

乙

(三)、程序與結果

1.二帶濾波器組

(1)源程序：

clear;

elf;

Fp=1700;Fr=2300;Fs=8000;

Wp=tan(pi*Fp/Fs);

Wr=tan(pi*Fr/Fs);

Wc=sqrt(Wp*Wr);

k=Wp/Wr;

kl=sqrt(sqrt(l-kA2));

q0=0.5*(1-kl)/(1+kl);

q=q0+2*q0A5+15*q0A9+150*q0A13;

N=ll;

N2=fix(N/4);

M=fix(N/2);

N1=M-N2;

forjj=1:M

a=0;

form=0:5

a=a+(-1)Am*qA(m*(m+1))*sin((z*m+l)*pi*jj/N);%Nis

odd,u=j

end

a

b=0;

form=l:5

b=b+(-1)Am*qA(mA2)*cos(2*m*pi*jj/N);

end

b

W(jj)=2*qA0.25*a/(l+2*b);

V(jj)=sqrt((l-k*W(jj)A2)*(1-W(jj)A2/k));

end

fori=l:N1

alpha(i)=2*V(2*i-l)/(l+W(2*i-l)A2);

end

fori=l:N2

beta(i)=2*V(2*i)/(1+W(2*i)A2);

end

fori=l:N1

a(i)=(1-alpha(i)*Wc+WcA2)/(1+alpha(i)*Wc+Wczs2);

end

fori=l:N2

b(i)=(1-beta(i)*Wc+WcA2)/(1+beta(i)*Wc+WcA2);

end

w=O:0.0001:0.5;

LP=zeros(size(w));HP=zeros(size(w));

forn=l:length(w)

z=exp(j*w(n)*2*pi);

Hl=l;

fori=l:N1

(a(i)+z^(-2))/(1+a(i)*zA(-2));

end

H2=l/z;

fori=l:N2

H2=H2*(b(i)+z"(-2))/(1+b(i)*zA(-2));

end

LP(n)=abs((H1+H2)/2);

HP(n)=abs((H1-H2)/2);

end

plot(w,LP,'r');

holdon;

xlabel(*digitalfrequency*);

ylabel(*amptitude');

(2)運行結果：見圖1

a

p

n

三

d

E

e

0

00.050.10.150.20.250.30.350.40.450.5

digitalfrequency

圖1二帶數字淀波器組

2.四帶濾波器組

(1)源程序：

clear;

elf;

Fp=1700;Fr=2300;Fs=8000;

Wp=tan(pi*Fp/Fs);

Wr=tan(pi*Fr/Fs);

Wc=sqrt(Wp*Wr);

k=Wp/Wr;

kl=sqrt(sqrt(l-kA2));

q0=0.5*(1-kl)/(1+kl);

q=qO+2*qOA5+15*qOA9+150*qOA13；

N=ll；

N2=fix(N/4);

M=fix(N/2);

N1=M-N2;

forjj=1:M

a=0;

form=0:5

a=a+(-1)AmxqA(m*(m+1))*sin((z*m+l)*pi*jj/N);%Nis

odd,u=j

end

b=0;

form=l:5

b=b+(-1)AmTqA(mA2)*cos(2*m*pi*jj/N);

end

W(jj)=2*qA0.25*a/(l+2*b);

V(jj)=sqrt((l-k*W(jj)A2)*(1-W(jj)A2/k));

End

fori=l:N1

alpha(i)=2*V(2*i-l)/(1+W(2*i-1)A2);

end

tori=l:N2

beta(i)=2*V(2*i)/(1+W(2*i)A2);

end

fori=l:N1

a(i)=(1-alpha(i)*Wc+WcA2)/(1+alpha(i)*Wc+WcA2);

end

fori=l:N2

A

b(i)=(1-beta(i)*Wc+Wc人2)/(1+beta(i)*Wc+Wc2);

end

w=0:0.0001:0.5;

LLP=zeroa(size(w));LHP=zero5(size(w));

HLP=zeros(size(w));HHP=zeros(size(w));

forn=l:length(w)

z=exp(j*w(n)*2*pi);

Hl=l;

fori=l:N1

H1=H1*(a(i)+zA(-2))/(l+a(i)*zA(-2));

end

H21=l;

fori=l:N1

H21=H21*(a(i)+zA(-4))/(l+a(i)*zA(-4));

end

H2=l/z;

fori=l:N2

H2=H2*(b(i)+zA(-2))/(1+b(i)*zA(-2));

end

H22=l/(z人2);

fori=l:N2

H22=H22*(bii)+zA(-4))/(l+b(i)*z^(-4));

end

LP=((H1+H2)/2);

HP=((H1-H2)/2);

LLP(n)=abs((H21+H22)/2*LP);

LHP(n)=abs((H21-H22)/2*LP);

HHP(n)=abs((H21+H22)/2*HP);

HLP(n)=abs((H21-H22)/2*HP);

end

plot(w,LLP,'bIw,LHP,'rIw,HLP,，kIw,HHP,1m')

holdon

xlabel(*digitalfrequency*);

ylabel(*amptitude1);

(2)運行結果：見圖2

9

8

o.7

o.

o,6

a

p

no.

三6

d

Eo.

eo.4

o.

o.'3

o.2

1

0

00.050.10.150.20.250.30.350.40.450.5

digitalfrequency

圖2四帶數字濾波器組

二、根據《現代數字信號處理》第四章提供的數據，試用如下方法估計其功

率譜，并畫出不同參數情況下的功率譜曲線：

1)Levison算法

2)Burg算法

3）ARMA模型法

4）MUSIC算法

1Levinson算法

Levinson算法用于求解Yule-Wa珠er方程，是一種按階次進行遞推的

算法，即首先以AR（0）和AR（1）模型參數作為初始條件，計算AR（2）

模型參數；然后根據這些參數計算AR（3）參數，等等，一直到計算出AR

（P）模型參數為止，需要的運算量數量級為〃2,其中p為AR模型的階數。

它利用了方程系數矩陣的對稱性和Toeplitz性質，是一種高效的算法。

Levinson算法的優點是計算簡單，步步監視均方誤差，其缺點是需要由觀

測數據計算自相關值.當觀測數據短時，誤差較大；當觀測數據長時，計算

量大；并且會產生譜峰漂移和譜線分裂。

算法步驟如下：

（1）由輸入數據估計自相關函數，一種漸近無偏估計（稱之為取樣自

相關函數）的公式為：

1N-1-卜

=—Zx*(n)x(m+n),\n\<N-\

N〃=o

（2）根據估計所得到的自相關函數，用下面的迭代公式估算AR模型參

數：

以=~0）+

f=l

2=Z4“R（A+j），4,0=°

/=0

哈=（1一|九+『）吠

=ak,i~/k+iak,k+\-i^'=1,2,…，Z

做+W+1——九+1

（3）對于AR（p）模型，按以上述遞推公式迭代計算直到左+1=〃時為

止。將算出的模型參數代入下式即可得到功率譜估計值：

2')=

1=1

Levinson算法的MATLAB源程序如下：

其中，參數X為輸入序列，p為AR模型的階數，函數調用形式為：Levinson

(X,p)0

functionS=Levinson(X,p)

N=length(X);

form=0:N-l

R(m+l)=sum(conj(X([1:N-m])).*X(fm+l:N]))/N;

end

a=-R(2)/R(l);

sigma2=(1-abs(a)A2)*R(1);

gamma=-a;

fork=2:p

sigma2(k)=R(1)+sum(a.*conj(R([2:k])));

D=R(k+l)+sum(a.*R([k:-1:2]));

gamma(k)=D/sigma2(k);

sigma2(k)=(1-abs(gamma(k))八2)*sigma2(k);

a=[a-gamma(k)*conj(fliplr(a)),-gamma(k)];

end

sigma2=real(signa2);

f=linspace(-0.5,0.5,512);

fork=l:512

S(k)=10*logl0(sigma2(end)/(abs(1+sum(a.*exp(-j*2*pi*f(

k)*[l：p])))A2));

end

plot(f,S);

title([*Levinson:*,int2str(p),*階']);

xlabel(1歸一化頻率1),ylabel「相對譜/dB，);

分別對應于10階、25階、40階、55階的Levinson算法的源程序和執

行結果如下所示：

elf

p=[10254055];

fork=l:4

subplot(2,2,k);

Levinson(xzp(k));

end

2O

1O

-1O

-29O

圖3lenxison算法

2Burg算法

Burg算法一方面希望利用已知數據段兩端以外的未知數據(它對這些

未知數據不作主觀臆測)，另一方面又總是設法保證使預測誤差濾波器是最

小相位的。Burg算法與自相關法和協方差法不同，它不直接計算AR參數,

而是先估計反射系數，然后利用Levinson遞推算法由反射系數求得AR參數，

估計反射系數時所使用的準則是前向和后向預測誤差功率估計的平均值最

小準則，在這里，預測誤差功率估計仍然用時間平均來代替集合平均。Burg

法估計反射系數的準則表示為：

N-\

￡二Z1e；（〃）25）2]=min

n=p

其優點是實現簡單，在一定程度上可以克服Levinson算法的譜峰漂移和譜

線分裂的缺點。

算法步驟如下：

(1)設輸入數據序列為x(〃),0W〃WN-1,對前后向預測誤差之和求偏

導，得反射系數

N-]

22叱|(〃)唱(〃-1)

y—k_____________________

/k一N-l22

Z(*(〃)l+肩(-1)|)

n=k

前后向預測誤差遞推公式:

回叫=(1-九丫、

(W)J〔一汽1人％5TL

ak,i=ak-\,i~Ykak-\,k-i^=1,2,3,...,匕。「\,0=1

(2)重復以上步驟直至七p,根據迭代得到的AR模型參數計算功率譜，

計算功率譜的公式同上面算法。

Burg算法的MATLAB源程序如下：

其中，參數X為輸入序列，p為AR模型的階數，函數調用形式為：Burg

(X,p)o

functionS=Burg[(X,p)

N=length(X);

ef=X;eb=X;

sigma2=sum(X*X')/N;

a=[];

fork=l:p

efp=ef(k+1:end);

ebp=eb(k:end-1);

gamma(k)=2*efp*ebp1/(efp*efp*+ebp*ebp1);

sigma2(k+1)=(1-abs(gamma(k))A2)*sigma2(k);

ef(k+1:end)=efp-gamma(k)*ebp;

eb(k+1：end)=ebp-gamma(k)**efp;

a=[a-gamma(k)*conj(fliplr(a)),-gamma(k)];

end

f=linspace(-0.5,0.5,512);

fork=l:512

S(k)=10*logl0(sigma2(end)/(abs(l+sum(a.*exp(-j*2*pi*f(

k)*[l:p])))A2));

end

plot(f,S);

title([*Burg:*)int2str(p),*階」)；

xlabel(，歸一化頻率,)，ylabel「相對譜/dB，);

分別對應于10階、25階、40階、55階的Burg算法的源程序和執行結

果如下所示：

elf

p=[10254055];

fork=l:4

subplot(2,2,k);

Burgl(x,p(k));

end

階

O

1

9

階

25

Burg:

BU

o一

3

40

o

2

2o

s

m

p

o

?

1、

^出

菽o

o

s

?哭

O

2

/

O

/

.1

0

40

_

-2

9

96

?5。

喟就喟0.5

4

o一

2

o’

m

m

40

p

p

20

痘0，、

器菽o

更2—

*

o

I

-20

_

-4

-40

5。

0

-0.5

歸化頻率0.5

率

化頻

歸一

算法

Burg

圖4

法

模型

RMA

3A

分

的高

可靠

得到

，能

g法時

Bur

采用

別是

法，特

計方

譜估

R模型

用A

當采

譜

好的

得良

能獲

型才