数字信号处理文档格式.docx

数字信号处理文档格式.docx

《数字信号处理文档格式.docx》由会员分享，可在线阅读，更多相关《数字信号处理文档格式.docx（25页珍藏版）》请在冰豆网上搜索。

rxx=xcorr（x）;

%自相关函数估计

fori=1:

100

forj=1:

mrxx（i,j）=rxx（1000-i+j）;

%维纳霍夫方程

end

end

xd=cos（sita）;

%产生维纳滤波中x方向上观测信号与期望信号的互相关矩阵

rxd=xcorr（x,xd）;

mrxd（i）=rxd（999+i）;

hoptx=inv（mrxx）*mrxd'

;

%由维纳-霍夫方程得到的x方向上的滤波器最优解（）

fx=conv（x,hoptx）;

%滤波后x方向上的输出

eminx=xd-fx（1:

1000）;

%x方向上最小均方误差

%产生维纳滤波中y方向上观测信号的自相关矩阵

ryy=xcorr（y）;

mryy（i,j）=ryy（1000-i+j）;

yd=sin（sita）;

%产生维纳滤波中y方向上观测信号与期望信号的互相关矩阵

ryd=xcorr（y,yd）;

mryd（i）=ryd（999+i）;

hopty=inv（mryy）*mryd'

%由维纳-霍夫方程得到的y方向上的滤波器最优解

fy=conv（y,hopty）;

%滤波后y方向上的输出

eminy=yd-fy（1:

%y方向上最小均方误差

figure;

subplot（2,2,1）

plot（xd）;

title（'

x方向期望信号'

）;

subplot（2,2,2）

plot（xnoise）;

x方向噪声信号'

subplot（2,2,3）

plot（x）;

x方向观测信号'

subplot（2,2,4）%输出x的滤波后信号

plot（fx）;

x滤波后信号'

plot（yd）;

y方向期望信号'

plot（ynoise）;

y方向噪声信号'

plot（y）;

y方向观测信号'

subplot（2,2,4）%输出y的滤波后信号

plot（fy）;

滤波后信号'

plot（xd,yd,'

k'

holdon;

plot（x,y,'

b:

'

plot（fx,fy,'

g'

最终结果'

3、运行结果展示：

4、遇到的问题和解决方法：

a.问：

在生成自相关矩阵时，循坏结构mrxx（i,j）=rxx（1000-i+j）中的数1000不能被任何数替代，这是为什么？

答：

因为在生成自相关函数时，生成2N-1个数（N=1000），为了生成的维纳霍夫方程

是对称的矩阵，生自相关函数的第1000个数正好对应维纳霍夫方程矩阵的第1

个数，即就是rxx（0）。

b.问：

X和Y方向上最小均方误差为eminx=xd-fx（1:

1000），此时的滤波器长度为1099，为什么当输入其他1000个长度时，画不出图像？

答：

因为xd期望信号是1:

1000的维度，为了使得与xd维度相匹配，故而要求fx一定是1:

1000，否则维度不匹配，无法进行运算。

C.问：

在最后一幅图中，绿色线表示滤波后的轨迹，在运动轨迹的中心有其运动轨迹的原因？

可能与滤波中的卷积运算有关，卷积使得信号长度增加，从而出现中心点，看起来有点多余。

二、（自适应滤波器设计）

自适应滤波器是通过最陡下降法和Widrow-HoffLMS算法来确定最佳权系数wj,以及通过LMS算法确定收敛系数μ，通过设定次数迭代，得到最佳滤波信号。

2、程序代码

clearall;

N=1500;

%采样点数

v=normrnd（0,6^0.5,1,N）;

%噪声信号

n=1:

N;

d=6*sin（0.03*n）;

%期望信号

dn=reshape（d,N,1）;

%从生成的期望信号中采出一个N行1列的列向量来形成dn矩阵

x=d+v;

%期望信号与噪声叠加后的观测信号1*N

X=reshape（x,N,1）;

%从实际信号中采样得到X矩阵对应的列向量

%X=x'

%W=zeros（N,1）;

%初始化权系数矩阵，将这个列向量赋0

y=zeros（N,1）;

%初始化输出信号矩阵，将列向量全赋0

M=200;

%滤波器的阶数

r=max（eig（X*X.'

））;

%输入信号相关矩阵的最大特征值

u=（1/r）;

%收敛因子0<

u<

1/r

en=zeros（M,1）;

%误差序列,en（k）表示第k次迭代时预期输出与实际输入的误差

ee=zeros（M,1）;

%均方误差，用以判别输出信号对期望信号的偏离程度，从而度量系统性能的好坏

W=zeros（M,N）;

%每一行代表一个加权参量,每一列代表-次迭代,初始为0

fork=M:

N

x=X（k:

-1:

k-M+1）;

%滤波器M个抽头的输入是指从X中按照倒序取出从k开始的M个样点M*1

y=W（:

k-1）.'

*x;

%滤波器的输出w为M*1

en（k）=dn（k）-y;

%第k次迭代的误差

W（:

k）=W（:

k-1）+2*u*en（k）*x;

%滤波器权值计算的迭代式最陡下降法的递推公式

ee（k）=en（k）*en（k）;

%求en的平方，以便于求均方误差

yn=ones（size（X））;

%

length（X）

yn（k）=W（:

end）.'

%将W中的最后一列用来计算

%图像输出

figure

（1）

holdon

plot（X,'

plot（dn,'

r'

grid;

观测信号（绿）和期望信号（红）'

figure

（2）

plot（yn,'

b'

输出信号（蓝）和期望信号（红）'

最优滤波信号'

figure（3）

subplot（211）

plot（ee）,grid;

均方误差'

subplot（212）

plot（W）;

权系数'

）

4、分析：

1、在第二幅图中，输出信号（蓝）的0-200内信号显示为0，这是因为滤波器的阶数为200，且求权值系数的循坏是从200-1000，故而前200信号显示为0.

2、yn=ones（size（X））输出最优滤波信号前面也可以加上inf，表示输出序列是指生成与X维数相同的一个矩阵，矩阵的每个值都是无穷大。

三、（功率谱）

经典谱估计法可以通过BT法得到理想功率谱，即就是先按照有限个观测数据估计自相关函数，再计算功率谱。

周期图法即直接对观测数据进行处理，计算功率谱，计算公式为：

closeall;

Fs=500;

%采样率

N=1024;

%观测数据

w=sqrt（0.36）+randn（1,N）;

%0均值,方差为0.36的白噪声,长度1024

x=[w

（1）zeros（1,N-1）];

%初始化x（n）,长度1024,x

（1）=w

（1）剩下的都为0

index=0:

round（N/2-1）;

fori=2:

x（i）=0.5*x（i-1）+w（i）;

%迭代产生观测数据x（n）

x1（i）=0.5*x（i-1）;

figure

plot（w）;

gridon;

wn'

xn'

%%理想功率谱

Pxx1=abs（fft（x1））;

Pxx1=10*log10（Pxx1（index+1））;

%化为db

plot（Pxx1）;

理想功率谱'

xlabel（'

频率'

ylabel（'

功率db'

%%周期图法

%1024个观测点

Pxx=abs（fft（x））.^2/N;

%周期图公式

%Pxx=10*log10（Pxx）;

Pxx=10*log10（Pxx（index+1））;

%f=（0:

length（Pxx）-1）*Fs/length（Pxx）;

%给出频率序列

plot（Pxx）;

周期图1024点'

%周期图256个观测点

x1=x（1:

4:

N）;

Pxx1=abs（fft（x1,1024））.^2/N;

%Pxx1=10*log10（Pxx1）;

周期图256点'

3、结果展示：

周期图法得到的功率谱估计，谱线的起伏较大，即估计所得的均方误差较大。

当N增加时，摆动的频率随N的变化而加快，而摆动的幅度变化并不是很大。

且N=1024时，谱的分辨率较N=256时大。

四、（序列的自功率谱和互功率谱）

利用循环迭代产生所需要的序列，然后先求得自相关函数，然后求得自动率谱。

clc;

X=sqrt

（1）+randn（1,N）;

%0均值,方差为1的白噪声,长度1024

Y=[X

（1）zeros（1,N-1）];

Z=[Y

（1）zeros（1,N-1）];

Fs=1000;

forn=2:

Y（n）=（X（n）+X（n-1）