数字信号处理.docx
- 文档编号:1905071
- 上传时间:2022-10-25
- 格式:DOCX
- 页数:25
- 大小:304.33KB
数字信号处理.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);%收敛因子0
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的平方,以便于求均方误差
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)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理