数字信号处理实验

实验一 自适应滤波器

一、实验目的

1、掌握功率谱估计方法

2、会用matlab 对功率谱进行仿真 二、实验原理

功率谱估计方法有很多种，一般分成两大类，一类是经典谱估计；另一类是现代谱估计。经典谱估计可以分成两种，一种是BT 法，另一种是周期法；BT 法是先估计自相关函数，然后将相关函数进行傅里叶变换得到功率谱函数。相应公式如下所示：

||1

*0

1

ˆ()()()(11)

ˆˆ()(12)

N m xx n jwn BT

xx m r

m x n x n m N P r

m e --=∞

-=-∞

=+-=-∑

∑

周期图法是采用功率谱的另一种定义，但与BT 法是等价的，相应的功率谱估计如下所示：

2

1

1ˆ

()()01

(13)N jw jwn xx n P e x n e n N N

--==≤≤--∑

其计算框图如下所示：

)

(jw xx e ∧

图1.1周期图法计算用功率谱框图

由于观测数据有限，所以周期图法估计分辨率低，估计误差大。针对经典谱估计的缺点，一般有三种改进方法：平均周期图法、窗函数法和修正的周期图平均法。

三、实验要求

信号是正弦波加正态零均值白噪声，信噪比为10dB，信号频率为2kHZ，取样频率为100kHZ。

四、实验程序与实验结果

（1）用周期图法进行谱估计

A、实验程序：

%用周期法进行谱估计

clear all;

N1=128;%数据长度

N2=256;

N3=512;

N4=1024;

f=2;%正弦波频率，单位为kHZ

fs=100;%抽样频率,单位为kHZ

n1=0:N1-1;

n2=0:N2-1;

n3=0:N3-1;

n4=0:N4-1;

a=sqrt(20);%由信噪比为10dB计算正弦信号的幅度

wn1=randn(1,N1);xn1=a*sin(2*pi*f*n1./fs)+wn1; Pxx1=10*log10(abs(fft(xn1).^2)/N1);%周期法求功率谱f1=((0:length(Pxx1)-1))/length(Pxx1);

wn2=randn(1,N2);xn2=a*sin(2*pi*f*n2./fs)+wn2; Pxx2=10*log10(abs(fft(xn2).^2)/N2);

f2=((0:length(Pxx2)-1))/length(Pxx2);

wn3=randn(1,N3);xn3=a*sin(2*pi*f*n3./fs)+wn3; Pxx3=10*log10(abs(fft(xn3).^2)/N3);

f3=((0:length(Pxx3)-1))/length(Pxx3);

wn4=randn(1,N4);xn4=a*sin(2*pi*f*n4./fs)+wn4; Pxx4=10*log10(abs(fft(xn4).^2)/N4);

f4=((0:length(Pxx4)-1))/length(Pxx4);

subplot(2,2,1);

plot(f1,Pxx1);xlabel('频率');ylabel('功率(dB)'); title('功率谱Pxx，N=128');

subplot(2,2,2);

plot(f2,Pxx2);xlabel('频率');ylabel('功率(dB)'); title('功率谱Pxx，N=256');

subplot(2,2,3);

plot(f3,Pxx3);xlabel('频率');ylabel('功率(dB)'); title('功率谱Pxx，N=512');

subplot(2,2,4);

plot(f4,Pxx4);xlabel('频率');ylabel('功率(dB)'); title('功率谱Pxx，N=1024');

B、实验仿真结果：

（2）采用汉明窗，分段长度L=32，用修正的周期图求平均法进行谱估计

A：实验程序:

clear all;

N=512;%数据长度

Ns=32;%分段长度

f1=2;%正弦波频率，单位为kHZ

fs=100;%抽样频率,单位为kHZ

n=0:N-1;

a=sqrt(20);%由信噪比为10dB计算正弦信号的幅度

wn=randn(1,N);

xn=a*sin(2*pi*f1*n./fs)+wn;

w=hamming(32)';%汉明窗

Pxx1=abs(fft(w.*xn(1:32),Ns).^2)/norm(w)^2;

Pxx2=abs(fft(w.*xn(33:64),Ns).^2)/norm(w)^2;

Pxx3=abs(fft(w.*xn(65:96),Ns).^2)/norm(w)^2;

Pxx4=abs(fft(w.*xn(97:128),Ns).^2)/norm(w)^2;

Pxx5=abs(fft(w.*xn(129:160),Ns).^2)/norm(w)^2;

Pxx6=abs(fft(w.*xn(161:192),Ns).^2)/norm(w)^2;

Pxx7=abs(fft(w.*xn(193:224),Ns).^2)/norm(w)^2;

Pxx8=abs(fft(w.*xn(225:256),Ns).^2)/norm(w)^2;

Pxx9=abs(fft(w.*xn(257:288),Ns).^2)/norm(w)^2;

Pxx10=abs(fft(w.*xn(289:320),Ns).^2)/norm(w)^2;

Pxx11=abs(fft(w.*xn(321:352),Ns).^2)/norm(w)^2;

Pxx12=abs(fft(w.*xn(353:384),Ns).^2)/norm(w)^2;

Pxx13=abs(fft(w.*xn(385:416),Ns).^2)/norm(w)^2;

Pxx14=abs(fft(w.*xn(417:448),Ns).^2)/norm(w)^2;

Pxx15=abs(fft(w.*xn(449:480),Ns).^2)/norm(w)^2;

Pxx16=abs(fft(w.*xn(481:512),Ns).^2)/norm(w)^2;

Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7+Pxx8+Pxx9 +Pxx10+Pxx11+Pxx12+Pxx13+Pxx14+Pxx15+Pxx16)/16);