假设一个点目标在x,y平面上绕单位圆做圆周运动,由于外界干扰,其运动轨迹发生了偏移。其中,x方向的干扰为均值为0,方差为0.05的高斯噪声;y方向干扰为均值为0,方差为0.06的高斯噪声。
1、
试设计一FIR维纳滤波器,确定最佳传递函数,并用该滤波器处理观测信号,得到其最佳估计。(注:自行设定误差判定阈值,根据阈值确定滤波器的阶数或传递函数的长度)。
4、 附件的实验报告中给出了解题思路,实现源程序、以及结果分析,
分别绘制了x方向和y方向的期望信号、噪声信号、观测信号、滤波后信号、最小均方误差信号的曲线图;
6、绘制了期望信号、观测信号和滤波后点目标的运动轨迹。
clear all
clc
k=8;
fs = 1000; %采样率
N = 1000; %采样点数
n = 0:N-1;
t = 0:1/fs:1-1/fs; %时间序列
signal=sin(2*pi*10*t);
noise=sin(2*pi*100*t); %前500点高斯分部白噪声,后500点均匀分布白噪声
xn= signal+k*noise; %构造的混合信号
figure(1)
plot(xn); grid on;
title('原始信号 ');
Fs=fft(xn,512); %将信号变换到频域
AFs=abs(Fs); %信号频域图的幅值
f=(0:255)*fs/512; %频率采样
figure(2)
plot(f,AFs(1:256)); %滤波前的信号频域图
grid on;
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
figure(3);
Rxx=xcorr(xn,xn); %得到混合信号的自相关函数
M=100; %维纳滤波器阶数
for i=1:M %得到混合信号的自相关矩阵
for j=1:M
rxx(i,j)=Rxx(abs(j-i)+N);
end
end
Rxy=xcorr(xn,signal); %得到混合信号和原信号的互相关函数
for i=1:M
rxy(i)=Rxy(i+N-1);
end %得到混合信号和原信号的互相关向量
h = inv(rxx)*rxy'; %得到所要涉及的wiener滤波器系数
Signal_Filter=filter(h,1, xn); %将输入信号通过维纳滤波器
plot(Signal_Filter);grid on;
title('维纳滤波后的信号');
Fs2=fft(Signal_Filter,512); %将信号变换到频域
AFs2=abs(Fs2); %信号频域图的幅值
f=(0:255)*fs/512; %频率采样
figure(4)
plot(f,AFs2(1:256)); %滤波前的信号频域图
grid on;
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
figure(5)
subplot(221)
plot(xn);grid on;
title('原始信号 ');
subplot(222)
plot(f,AFs(1:256)); %滤波前的信号频域图
grid on;
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图');
subplot(223)
plot(Signal_Filter);grid on;
title('维纳滤波后的信号');
subplot(224)
plot(f,AFs2(1:256)); %滤波前的信号频域图
xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');
grid on;
%实例一程序_维纳滤波的计算机实现
%初步处理,并接受输入数据,包括信号样本个数L和滤波器阶数N
clear all
close all
L=input('L=');
N=input('N=');
a=0.95;
%定义w,v,u
w=sqrt(3*(1-a^2))*(2*rand(1,L)-1);
v=sqrt(3)*(2*rand(1,L)-1);
u=ones(1,L);
%获得原信号s(n)和带噪声信号x(n)
s(1)=1;
for i=2:L
s(i)=a*s(i-1)+w(i);
end
for i=1:L
x(i)=s(i)+v(i);
end
%绘图比较原信号s(n)和带噪声信号x(n)
figure
k=(L-99):L;
plot(k,s(k),'r',k,x(k),'b');
legend('s(n)','x(n)',0);
title('comparation between s(n) and x(n)');
xlabel('n');ylabel('Input');
%计算信号x(n)的N阶自相关矩阵Rxx,x(n)与s(n)的互相关函数向量rxs,并获得估计FIR滤波器系数h1
phixx=xcorr(x,x);
for i=1:N
for j=1:N
Rxx(i,j)=phixx(i-j+L);
end
end
phixs=xcorr(x,s);
for i=1:N
rxs(i)=phixs(i+L);
end
h1=(inv(Rxx))*rxs';
%获得理想FIR滤波器系数h1
for i=1:N
h(i)=0.238*0.724^i*u(i);
end
%实例一程序
%绘图比较估计滤波器与理想滤波器
figure
k=1:N;
plot(k,h(k),'r',k,h1(k),'b');
title('Ideal h(n) & Calculated h(n)');
legend('Ideal h(n)',' Calculated h(n)');
xlabel('n');ylabel('h(n)');
%计算并绘图比较理想输出与实际输出
S=conv(h,v);
SI(1)=S(1);
for i=2:L
SI(i)=0.724*SI(i-1)+0.238*x(i);
end
figure
k=(L-99):L;
plot(k,s(k),'r',k,SI(k),'b');
title('s(n) VS. SI(n)');
legend('s(n)','SI(n)',0);
xlabel('n');ylabel('Ideal Output');
SR=conv(h1,x);
figure
k=(L-99):L;
plot(k,s(k),'r',k,SR(k),'b');
title('s(n)VS. SR(n)');
legend('s(n)','SR(n)',0);
xlabel('n');ylabel('Actual Output');
EX2=0;
EI2=0;
ER2=0;
%计算并输出所获得信号与原信号的均方误差,理想维纳滤波和估计维纳滤波均方误差
for i=1:L;
EX2=1/L*((x(i)-s(i))^2)+EX2;
end;
for i=1:L;
EI2=1/L*((SI(i)-s(i))^2)+EI2;
end;
for i=1:L;
ER2=1/L*((SR(i)-s(i))^2)+ER2;
end;
EX2
EI2
ER2