数字信号处理.docx
《数字信号处理.docx》由会员分享，可在线阅读，更多相关《数字信号处理.docx（25页珍藏版）》请在冰豆网上搜索。
数字信号处理
现代数字信号处理
实验报告
一、（设计FIR维纳滤波器）
1、解题思路：
由于题目中要求分别求X和Y方向上的信号的观测、噪声、输出等各种图像，并且在X和Y方向上的噪声也不同，故而就将信号分为X和Y轴。
由维纳霍夫方程的矩阵形式为了求hopt（=Rxx-1Rxd）,我们必须先求出观测信号的自相关矩阵RXX互相关矩阵Rxd，在定义时要注意选择相关函数的无偏性。
再由最小均方误差定义求出ex。
y方向上的计算与x方向相一致。
最后使x方向和y方向上的信号分别通过最优滤波，再进行合成，从而最终得到最优滤波后的观测信号。
2、程序代码：
clear;
clf;
sita=0:
pi/499.5:
2*pi;
xnoise=sqrt（0.03）*randn（1,1000）;%产生x轴方向噪声
ynoise=sqrt（0.02）*randn（1,1000）;%产生y轴方向噪声
x=cos（sita）+xnoise;%产生x轴方向观测信号
y=sin（sita）+ynoise;%产生y轴方向观测信号
%产生维纳滤波中x方向上观测信号的自相关矩阵
rxx=xcorr（x）;%自相关函数估计
fori=1:
100
forj=1:
100
mrxx（i,j）=rxx（1000-i+j）;%维纳霍夫方程
end
end
xd=cos（sita）;
%产生维纳滤波中x方向上观测信号与期望信号的互相关矩阵
rxd=xcorr（x,xd）;
fori=1:
100
mrxd（i）=rxd（999+i）;
end
hoptx=inv（mrxx）*mrxd';%由维纳-霍夫方程得到的x方向上的滤波器最优解（）
fx=conv（x,hoptx）;%滤波后x方向上的输出
eminx=xd-fx（1:
1000）;%x方向上最小均方误差
%产生维纳滤波中y方向上观测信号的自相关矩阵
ryy=xcorr（y）;
fori=1:
100
forj=1:
100
mryy（i,j）=ryy（1000-i+j）;
end
end
yd=sin（sita）;
%产生维纳滤波中y方向上观测信号与期望信号的互相关矩阵
ryd=xcorr（y,yd）;
fori=1:
100
mryd（i）=ryd（999+i）;
end
hopty=inv（mryy）*mryd';%由维纳-霍夫方程得到的y方向上的滤波器最优解
fy=conv（y,hopty）;%滤波后y方向上的输出
eminy=yd-fy（1:
1000）;%y方向上最小均方误差
figure;
subplot（2,2,1）
plot（xd）;
title（'x方向期望信号'）;
subplot（2,2,2）
plot（xnoise）;
title（'x方向噪声信号'）;
subplot（2,2,3）
plot（x）;
title（'x方向观测信号'）;
subplot（2,2,4）%输出x的滤波后信号
plot（fx）;
title（'x滤波后信号'）;
figure;
subplot（2,2,1）
plot（yd）;
title（'y方向期望信号'）;
subplot（2,2,2）
plot（ynoise）;
title（'y方向噪声信号'）;
subplot（2,2,3）
plot（y）;
title（'y方向观测信号'）;
subplot（2,2,4）%输出y的滤波后信号
plot（fy）;
title（'滤波后信号'）;
figure;
plot（xd,yd,'k'）;
holdon;
plot（x,y,'b:
'）;
holdon;
plot（fx,fy,'g'）;
title（'最终结果'）;
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.问：
在最后一幅图中，绿色线表示滤波后的轨迹，在运动轨迹的中心有其运动轨迹的原因？
答：
可能与滤波中的卷积运算有关，卷积使得信号长度增加，从而出现中心点，看起来有点多余。
二、（自适应滤波器设计）
1、解题思路：
自适应滤波器是通过最陡下降法和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）;%收敛因子0en=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的平方，以便于求均方误差
end
yn=ones（size（X））;%
fork=M:
length（X）
x=X（k:
-1:
k-M+1）;
yn（k）=W（:
end）.'*x;%将W中的最后一列用来计算
end
%图像输出
figure
（1）
holdon
plot（X,'g'）;
plot（dn,'r'）;grid;
title（'观测信号（绿）和期望信号（红）'）;
figure
（2）
holdon
plot（yn,'b'）;
plot（dn,'r'）;grid;
title（'输出信号（蓝）和期望信号（红）'）;
title（'最优滤波信号'）;
figure（3）
subplot（211）
plot（ee）,grid;
title（'均方误差'）;
subplot（212）
plot（W）;grid;
title（'权系数'）
3、运行结果展示：
4、分析：
1、在第二幅图中，输出信号（蓝）的0-200内信号显示为0，这是因为滤波器的阶数为200，且求权值系数的循坏是从200-1000，故而前200信号显示为0.
2、yn=ones（size（X））输出最优滤波信号前面也可以加上inf，表示输出序列是指生成与X维数相同的一个矩阵，矩阵的每个值都是无穷大。
三、（功率谱）
1、解题思路：
经典谱估计法可以通过BT法得到理想功率谱，即就是先按照有限个观测数据估计自相关函数，再计算功率谱。
周期图法即直接对观测数据进行处理，计算功率谱，计算公式为：
2、程序代码：
clear;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:
N
x（i）=0.5*x（i-1）+w（i）;%迭代产生观测数据x（n）
end
fori=2:
N
x1（i）=0.5*x（i-1）;%迭代产生观测数据x（n）
end
figure
subplot（211）
plot（w）;
gridon;
title（'wn'）;
subplot（212）
plot（x）;
title（'xn'）;
gridon;
%%理想功率谱
Pxx1=abs（fft（x1））;
Pxx1=10*log10（Pxx1（index+1））;%化为db
figure;
plot（Pxx1）;
gridon;
title（'理想功率谱'）;
xlabel（'频率'）;
ylabel（'功率db'）;
%%周期图法
%1024个观测点
Pxx=abs（fft（x））.^2/N;%周期图公式
%Pxx=10*log10（Pxx）;
Pxx=10*log10（Pxx（index+1））;%化为db
%f=（0:
length（Pxx）-1）*Fs/length（Pxx）;%给出频率序列
figure;
plot（Pxx）;
gridon;
title（'周期图1024点'）;
xlabel（'频率'）;
ylabel（'功率db'）;
%周期图256个观测点
x1=x（1:
4:
N）;
Pxx1=abs（fft（x1,1024））.^2/N;%周期图公式
%Pxx1=10*log10（Pxx1）;
Pxx1=10*log10（Pxx1（index+1））;%化为db
figure;
plot（Pxx1）;
gridon;
title（'周期图256点'）;
xlabel（'频率'）;
ylabel（'功率db'）;
3、结果展示：
4、分析：
周期图法得到的功率谱估计，谱线的起伏较大，即估计所得的均方误差较大。
当N增加时，摆动的频率随N的变化而加快，而摆动的幅度变化并不是很大。
且N=1024时，谱的分辨率较N=256时大。
四、（序列的自功率谱和互功率谱）
1、解题思路：
利用循环迭代产生所需要的序列，然后先求得自相关函数，然后求得自动率谱。
2、程序代码：
clc;
clear;
closeall;
N=1024;
X=sqrt
（1）+randn（1,N）;%0均值,方差为1的白噪声,长度1024
Y=[X
（1）zeros（1,N-1）];%初始化x（n）,长度1024,x
（1）=w
（1）剩下的都为0
Z=[Y
（1）zeros（1,N-1）];
Fs=1000;
index=0:
round（N/2-1）;
forn=2:
N
Y（n）=（X（n）+X（n-1）
升级会员 