現(xiàn)代信號(hào)處理大作業(yè)

現(xiàn)代信號(hào)處理大型作業(yè)

一.試用奇階互補(bǔ)法設(shè)計(jì)兩帶濾波器組(高、低通互補(bǔ))，進(jìn)而實(shí)現(xiàn)四帶濾

波器組；并畫出其頻響。濾波器設(shè)計(jì)參數(shù)為：Fp=1.7KHz,F「=23KHz,Fs

=8KHz,加侖70dB0

(一)、分析

與通常的濾波器相比，互補(bǔ)濾波器具有優(yōu)良的結(jié)構(gòu)特性和結(jié)構(gòu)特性，具

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

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

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

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

=1,0<W<2K

則稱這組濾波器為功率互補(bǔ)濾波器，同時(shí)互補(bǔ)濾波器還應(yīng)該滿足：

XHk(z)=A(z)其中A(z)為全通函數(shù)，適當(dāng)?shù)倪x擇全通函數(shù)，可以

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

(二)、設(shè)計(jì)步驟

(1)對(duì)Fp、Fr進(jìn)行預(yù)畸

V\Fp

=爾（

Q=應(yīng)（

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

鏡像濾波器

(3)計(jì)算相關(guān)系數(shù)

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

,Qp

k="一人

ii-VF

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

N=lg(42/i6)/lgq；

i；(N活數(shù))

〃[，T(N為禺?dāng)?shù))

2公￡(_1尸產(chǎn)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)互補(bǔ)鏡像濾波器的數(shù)字實(shí)現(xiàn)

2-a,

A=

2+a.

2-A

B

i=2+0

口4+z々

H.(Z)=1,…M

1J1+AZ-2

pI7-2

區(qū)“心口七毛八卜以

j1TUj乙

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

乙

(三)、程序與結(jié)果

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)運(yùn)行結(jié)果：見圖1

a

p

n

三

d

E

e

0

00.050.10.150.20.250.30.350.40.450.5

digitalfrequency

圖1二帶數(shù)字淀波器組

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)運(yùn)行結(jié)果：見圖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四帶數(shù)字濾波器組

二、根據(jù)《現(xiàn)代數(shù)字信號(hào)處理》第四章提供的數(shù)據(jù)，試用如下方法估計(jì)其功

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

1)Levison算法

2)Burg算法

3）ARMA模型法

4）MUSIC算法

1Levinson算法

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

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

模型參數(shù)；然后根據(jù)這些參數(shù)計(jì)算AR（3）參數(shù)，等等，一直到計(jì)算出AR

（P）模型參數(shù)為止，需要的運(yùn)算量數(shù)量級(jí)為〃2,其中p為AR模型的階數(shù)。

它利用了方程系數(shù)矩陣的對(duì)稱性和Toeplitz性質(zhì)，是一種高效的算法。

Levinson算法的優(yōu)點(diǎn)是計(jì)算簡(jiǎn)單，步步監(jiān)視均方誤差，其缺點(diǎn)是需要由觀

測(cè)數(shù)據(jù)計(jì)算自相關(guān)值.當(dāng)觀測(cè)數(shù)據(jù)短時(shí)，誤差較大；當(dāng)觀測(cè)數(shù)據(jù)長(zhǎng)時(shí)，計(jì)算

量大；并且會(huì)產(chǎn)生譜峰漂移和譜線分裂。

算法步驟如下：

（1）由輸入數(shù)據(jù)估計(jì)自相關(guān)函數(shù)，一種漸近無(wú)偏估計(jì)（稱之為取樣自

相關(guān)函數(shù)）的公式為：

1N-1-卜

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

N〃=o

（2）根據(jù)估計(jì)所得到的自相關(guān)函數(shù)，用下面的迭代公式估算AR模型參

數(shù)：

以=~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）對(duì)于AR（p）模型，按以上述遞推公式迭代計(jì)算直到左+1=〃時(shí)為

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

2')=

1=1

Levinson算法的MATLAB源程序如下：

其中，參數(shù)X為輸入序列，p為AR模型的階數(shù)，函數(shù)調(diào)用形式為：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「相對(duì)譜/dB，);

分別對(duì)應(yīng)于10階、25階、40階、55階的Levinson算法的源程序和執(zhí)

行結(jié)果如下所示：

elf

p=[10254055];

fork=l:4

subplot(2,2,k);

Levinson(xzp(k));

end

2O

1O

-1O

-29O

圖3lenxison算法

2Burg算法

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

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

小相位的。Burg算法與自相關(guān)法和協(xié)方差法不同，它不直接計(jì)算AR參數(shù),

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

估計(jì)反射系數(shù)時(shí)所使用的準(zhǔn)則是前向和后向預(yù)測(cè)誤差功率估計(jì)的平均值最

小準(zhǔn)則，在這里，預(yù)測(cè)誤差功率估計(jì)仍然用時(shí)間平均來(lái)代替集合平均。Burg

法估計(jì)反射系數(shù)的準(zhǔn)則表示為：

N-\

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

n=p

其優(yōu)點(diǎn)是實(shí)現(xiàn)簡(jiǎn)單，在一定程度上可以克服Levinson算法的譜峰漂移和譜

線分裂的缺點(diǎn)。

算法步驟如下：

(1)設(shè)輸入數(shù)據(jù)序列為x(〃),0W〃WN-1,對(duì)前后向預(yù)測(cè)誤差之和求偏

導(dǎo)，得反射系數(shù)

N-]

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

y—k_____________________

/k一N-l22

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

n=k

前后向預(yù)測(cè)誤差遞推公式:

回叫=(1-九丫、

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

ak,i=ak-\,i~Ykak-\,k-i^=1,2,3,...,匕?！竆,0=1

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

計(jì)算功率譜的公式同上面算法。

Burg算法的MATLAB源程序如下：

其中，參數(shù)X為輸入序列，p為AR模型的階數(shù)，函數(shù)調(diào)用形式為：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「相對(duì)譜/dB，);

分別對(duì)應(yīng)于10階、25階、40階、55階的Burg算法的源程序和執(zhí)行結(jié)

果如下所示：

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法時(shí)

Bur

采用

別是

法，特

計(jì)方

譜估

R模型

用A

當(dāng)采

譜

好的

得良

能獲

型才