Kalman滤波和MATLAB实现
卡尔曼滤波器（Kalman Filter）是一个最优化自回归数据处理算法（optimal recursive data processing algorithm）。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的。他的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等
简单来说，卡尔曼滤波器是一个“optimal recursive data processing algorithm（最优化自回归数据处理算法）”。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的。他的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。
为了可以更加容易的理解卡尔曼滤波器，这里应用形象的描述方法讲解，不像参考书那样罗列一大堆的数学公式和数学符号。但是，他的5条公式是其核心内容。结合现代的计算机，其实卡尔曼的程序相当的简单，只要你理解了他的那5条公式。
假设我们要研究的对象是一个房间的温度。根据你的经验判断，这个房间的温度是恒定的，也就是下一分钟的温度等于现在这一分钟的温度。假设你对你的经验不是100%的相信，可能会有上下偏差几度。我们把这些偏差看成是高斯白噪声，也就是这些偏差跟前后时间是没有关系的而且符合高斯分配。另外，我们在房间里放一个温度计，但是这个温度计也是不准确的，测量值会比实际值有偏差。我们也把这些偏差看成是高斯白噪声。
好了，现在对于某一分钟我们有两个有关该房间的温度值：你根据经验的预测值（系统的预测值）和温度计的值（测量值）。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。
假如我们要估算k时刻的实际温度值。首先你要根据k-1时刻的温度值，来预测k时刻的温度。因为你相信温度是恒定的，所以你会得到k时刻的温度预测值是跟k-1时刻一样的，假定是23度，同时该值的高斯噪声的偏差是5度（5是这样得到的：如果k-1时刻估算出的最优温度值的偏差是3，你对自己预测的不确定度是4度，他们平方相加再开方，就是5）。然后，你从温度计那里得到了k时刻的温度值，假设是25度，同时该值的偏差是4度。
由于我们用于估算k时刻的实际温度有两个温度值，分别是23度和25度。究竟实际温度是多少呢？相信自己还是相信温度计呢？究竟相信谁多一点，我们可以用他们的协方差（covariance）来判断。因为Kg^2=5^2/(5^2+4^2)，所以Kg=0.78，我们可以估算出k时刻的实际温度值是：23+0.78*(25-23)=24.56度。可以看出，因为温度计的covariance比较小（比较相信温度计），所以估算出的最优温度值偏向温度计的值。
现在我们已经得到k时刻的最优温度值了，下一步就是要进入k+1时刻，进行新的最优估算。到现在为止，好像还没看到什么自回归的东西出现。对了，在进入k+1时刻之前，我们还要算出k时刻那个最优值（24.56度）的偏差。算法如下：((1-Kg)*5^2)^0.5=2.35。这里的5就是上面的k时刻你预测的那个23度温度值的偏差，得出的2.35就是进入k+1时刻以后k时刻估算出的最优温度值的偏差（对应于上面的3）。
就是这样，卡尔曼滤波器就不断的把covariance递归，从而估算出最优的温度值。他运行的很快，而且它只保留了上一时刻的covariance。上面的Kg，就是卡尔曼增益（Kalman Gain）。他可以随不同的时刻而改变他自己的值。
Dr Kalman的卡尔曼滤波器。涉及一些基本的概念知识，包括概率（Probability），随机变量（Random Variable），高斯或正态分配（Gaussian Distribution）等。
首先，要引入一个离散控制过程的系统。该系统可用一个线性随机微分方程来描述：
X(k)=A X(k-1)+B U(k)+W(k)
再加上系统的测量值：
Z(k)=H X(k)+V(k)
上两式子中，X(k)是k时刻的系统状态，U(k)是k时刻对系统的控制量。A和B是系统参数，对于多模型系统，他们为矩阵。Z(k)是k时刻的测量值，H是测量系统的参数，对于多测量系统，H为矩阵。W(k)和V(k)分别表示过程和测量的噪声。他们被假设成高斯白噪声，他们的协方差（covariance）分别是Q，R（这里我们假设他们不随系统状态变化而变化）。
对于满足上面的条件(线性随机微分系统，过程和测量都是高斯白噪声)，卡尔曼滤波器是最优的信息处理器。下面我们来用他们结合他们的covariances来估算系统的最优化输出（类似上一节那个温度的例子）。
首先我们要利用系统的过程模型，来预测下一状态的系统。假设现在的系统状态是k，根据系统的模型，可以基于系统的上一状态而预测出现在状态：
X(k|k-1)=A X(k-1|k-1)+B U(k) (1)
式(1)中，X(k|k-1)是利用上一状态预测的结果，X(k-1|k-1)是上一状态最优的结果，U(k)为现在状态的控制量，如果没有控制量，它可以为0。
到现在为止，我们的系统结果已经更新了，可是，对应于X(k|k-1)的covariance还没更新。我们用P表示covariance：
P(k|k-1)=A P(k-1|k-1)A’+Q (2)
式(2)中，P(k|k-1)是X(k|k-1)对应的covariance，P(k-1|k-1)是X(k-1|k-1)对应的covariance，A’表示A的转置矩阵，Q是系统过程的covariance。式子1，2就是卡尔曼滤波器5个公式当中的前两个，也就是对系统的预测。
现在我们有了现在状态的预测结果，然后我们再收集现在状态的测量值。结合预测值和
测量值，我们可以得到现在状态(k)的最优化估算值X(k|k)：
X(k|k)=X(k|k-1)+Kg(k)(Z(k)-H X(k|k-1)) (3)
其中Kg为卡尔曼增益(Kalman Gain)：
Kg(k)=P(k|k-1)H’/(H P(k|k-1)H’+R) (4)
到现在为止，我们已经得到了k状态下最优的估算值X(k|k)。但是为了要另卡尔曼滤波器不断的运行下去直到系统过程结束，我们还要更新k状态下X(k|k)的covariance：
P(k|k)=（I-Kg(k)H）P(k|k-1) (5)
其中I为1的矩阵，对于单模型单测量，I=1。当系统进入k+1状态时，P(k|k)就是式子(2)的P(k-1|k-1)。这样，算法就可以自回归的运算下去。
卡尔曼滤波器的原理基本描述了，式子1，2，3，4和5就是他的5个基本公式。根据这5个公式，可以很容易的实现计算机的程序。
这里举一个简单的例子来说明卡尔曼滤波器的工作过程。把房间看成一个系统，然后对这个系统建模。房间的温度是跟前一时刻的温度相同的，所以A=1。没有控制量，所以U(k)=0。因此得出：
X(k|k-1)=X(k-1|k-1) (6)
式子（2）可以改成：
P(k|k-1)=P(k-1|k-1)+Q (7)
因为测量的值是温度计的，跟温度直接对应，所以H=1。式子3，4，5可以改成以下：
X(k|k)=X(k|k-1)+Kg(k)(Z(k)-X(k|k-1)) (8)
Kg(k)=P(k|k-1)/(P(k|k-1)+R) (9)
P(k|k)=（1-Kg(k)）P(k|k-1) (10)
现在模拟一组测量值作为输入。假设房间的真实温度为25度，模拟200个测量值，这些测量值的平均值为25度，但是加入了标准偏差为几度的高斯白噪声。
为了令卡尔曼滤波器开始工作，需要告诉卡尔曼两个零时刻的初始值，是X(0|0)和P(0|0)。他们的值不用太在意，随便给一个就可以了，因为随着卡尔曼的工作，X会逐渐的收敛。但是对于P，一般不要取0，因为这样可能会令卡尔曼完全相信你给定的X(0|0)是系统最优的，从而使算法不能收敛。选取X(0|0)=1度，P(0|0)=10。
该系统的真实温度为25度，图中用黑线表示。图中红线是卡尔曼滤波器输出的最优化结果（该结果在算法中设置了Q=1e-6，R=1e-1）。
附matlab下面的kalman滤波程序：
clear
N=200;
w(1)=0;
w=randn(1,N)
x(1)=0;
a=1;
for k=2:N;
x(k)=a*x(k-1)+w(k-1);
end
V=randn(1,N);
q1=std(V);
Rvv=q1.^2;
q2=std(x);
Rxx=q2.^2;
q3=std(w);
Rww=q3.^2;
c=0.2;
Y=c*x+V;
p(1)=0;
s(1)=0;
for t=2:N;
p1(t)=a.^2*p(t-1)+Rww;
b(t)=c*p1(t)/(c.^2*p1(t)+Rvv);
s(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1));
p(t)=p1(t)-c*b(t)*p1(t);
end
t=1:N;
plot(t,s,'r',t,Y,'g',t,x,'b');
function[x,V,VV,loglik]=kalman_filter(y,A,C,Q,R,init_x,init_V,varargin) %Kalman filter.
%[x,V,VV,loglik]=kalman_filter(y,A,C,Q,R,init_x,init_V,...)
%
%INPUTS:
%y(:,t)-the observation at time t
%A-the system matrix
%C-the observation matrix
%Q-the system covariance
%R-the observation covariance
%init_x-the initial state(column)vector
%init_V-the initial state covariance
%
%OPTIONAL INPUTS(string/value pairs[default in brackets])
%'model'-model(t)=m means use params from model m at time t[ones(1,T)] %In this case,all the above matrices take an additional final dimension,
%i.e.,A(:,:,m),C(:,:,m),Q(:,:,m),R(:,:,m).
%However,init_x and init_V are independent of model(1).
%'u'-u(:,t)the control signal at time t[[]]
%'B'-B(:,:,m)the input regression matrix for model m
%
%OUTPUTS(where X is the hidden state being estimated)
%x(:,t)=E[X(:,t)|y(:,1:t)]
%V(:,:,t)=Cov[X(:,t)|y(:,1:t)]
%VV(:,:,t)=Cov[X(:,t),X(:,t-1)|y(:,1:t)]t>=2
%loglik=sum{t=1}^T log P(y(:,t))
%
%If an input signal is specified,we also condition on it:
%e.g.,x(:,t)=E[X(:,t)|y(:,1:t),u(:,1:t)]
%If a model sequence is specified,we also condition on it:
%e.g.,x(:,t)=E[X(:,t)|y(:,1:t),u(:,1:t),m(1:t)]
[os T]=size(y);
ss=size(A,1);%size of state space
%set default params
model=ones(1,T);
u=[];
B=[];
ndx=[];
args=varargin;
nargs=length(args);
for i=1:2:nargs
switch args
case'model',model=args{i+1};
case'u',u=args{i+1};
case'B',B=args{i+1};
case'ndx',ndx=args{i+1};
otherwise,error(['unrecognized argument'args]) end
end
x=zeros(ss,T);
V=zeros(ss,ss,T);
VV=zeros(ss,ss,T);
loglik=0;
for t=1:T
m=model(t);
if t==1
%prevx=init_x(:,m);
%prevV=init_V(:,:,m);
prevx=init_x;
prevV=init_V;
initial=1;
else
prevx=x(:,t-1);
prevV=V(:,:,t-1);
initial=0;
end
if isempty(u)
[x(:,t),V(:,:,t),LL,VV(:,:,t)]=...
kalman_update(A(:,:,m),C(:,:,m),Q(:,:,m),R(:,:,m),y(:,t),prevx,prevV,'initial',initial); else
if isempty(ndx)
[x(:,t),V(:,:,t),LL,VV(:,:,t)]=...
kalman_update(A(:,:,m),C(:,:,m),Q(:,:,m),R(:,:,m),y(:,t),prevx,prevV,...
'initial',initial,'u',u(:,t),'B',B(:,:,m));
else
i=ndx;
%copy over all elements;only some will get updated
x(:,t)=prevx;
prevP=inv(prevV);
prevPsmall=prevP(i,i);
prevVsmall=inv(prevPsmall);
[x(i,t),smallV,LL,VV(i,i,t)]=...
kalman_update(A(i,i,m),C(:,i,m),Q(i,i,m),R(:,:,m),y(:,t),prevx(i),prevVsmall,...
'initial',initial,'u',u(:,t),'B',B(i,:,m));
smallP=inv(smallV);
prevP(i,i)=smallP;
V(:,:,t)=inv(prevP);
end
end
loglik=loglik+LL;
end
最佳线性滤波理论起源于40年代美国科学家Wiener和前苏联科学家Ｋолмогоров等人的研究工作，后人统称为维纳滤波理论。从理论上说，维纳滤波的最大缺点是必须用到无限过去的数据，不适用于实时处理。为了克服这一缺点，60年代Kalman把状态空间模型引入滤波理论，并导出了一套递推估计算法，后人称之为卡尔曼滤波理论。卡尔曼滤波是以最小均方误差为估计的最佳准则，来寻求一套递推估计的算法，其基本思想是：采用信号与噪声的状态空间模型，利用前一时刻地估计值和现时刻的观测值来更新对状态变量的估计，求出现时刻的估计值。它适合于实时处理和计算机运算。
现设线性时变系统的离散状态方程和观测方程为：
X(k)=F(k,k-1)·X(k-1)+T(k,k-1)·U(k-1)
Y(k)=H(k)·X(k)+N(k)
其中
X(k)和Y(k)分别是k时刻的状态矢量和观测矢量
F(k,k-1)为状态转移矩阵
U(k)为k时刻动态噪声
T(k,k-1)为系统控制矩阵
H(k)为k时刻观测矩阵
N(k)为k时刻观测噪声
则卡尔曼滤波的算法流程为：
预估计X(k)^=F(k,k-1)·X(k-1)
计算预估计协方差矩阵
C(k)^=F(k,k-1)×C(k)×F(k,k-1)'+T(k,k-1)×Q(k)×T(k,k-1)'
Q(k)=U(k)×U(k)'
计算卡尔曼增益矩阵
K(k)=C(k)^×H(k)'×[H(k)×C(k)^×H(k)'+R(k)]^(-1)
R(k)=N(k)×N(k)'
更新估计
X(k)~=X(k)^+K(k)×[Y(k)-H(k)×X(k)^]
计算更新后估计协防差矩阵
C(k)~=[I-K(k)×H(k)]×C(k)^×[I-K(k)×H(k)]'+K(k)×R(k)×K(k)'
X(k+1)=X(k)~
C(k+1)=C(k)~
重复以上步骤
Ｍatlab实现代码
******************************************************************************* **************************************************
%%%%Constant Velocity Model Kalman Filter Simulation%%%%
%==================================================================== ======
clear all;close all;clc;
%%Initial condition
ts=1;%Sampling time
t=[0:ts:100];
T=length(t);
%%Initial state
x=[040020]';
x_hat=[0000]';
%%Process noise covariance
q=5
Q=q*eye(2);
%%Measurement noise covariance
r=5
R=r*eye(2);
%%Process and measurement noise
w=sqrt(Q)*randn(2,T);%Process noise
v=sqrt(R)*randn(2,T);%Measurement noise
%%Estimate error covariance initialization
p=5;
P(:,:,1)=p*eye(4);
%==================================================================== ======
%%Continuous-time state space model
%{
x_dot(t)=Ax(t)+Bu(t)
z(t)=Cx(t)+Dn(t)
%}
A=[0100;
0000;
0001;
0000];
B=[00;
10;
00;
01];
C=[1000;
0010];
D=[10;
01];
%%Discrete-time state space model
%{
x(k+1)=Fx(k)+Gw(k)
z(k)=Hx(k)+Iv(k)
Continuous to discrete form by zoh
%}
sysc=ss(A,B,C,D);
sysd=c2d(sysc,ts,'zoh');
[F G H I]=ssdata(sysd);
%%Practice state of target
for i=1:T-1
x(:,i+1)=F*x(:,i);
end
x=x+G*w;%State variable with noise
z=H*x+I*v;%Measurement value with noise
%==================================================================== ======
%%%Kalman Filter
for i=1:T-1
%%Prediction phase
x_hat(:,i+1)=F*x_hat(:,i);
%State estimate predict
P(:,:,i+1)=F*P(:,:,i)*F'+G*Q*G';
%Tracking error covariance predict
P_predicted(:,:,i+1)=P(:,:,i+1);
%%Kalman gain
K=P(:,:,i+1)*H'*inv(H*P(:,:,i+1)*H'+R);
%%Updata step
x_hat(:,i+1)=x_hat(:,i+1)+K*(z(:,i+1)-H*x_hat(:,i+1));
%State estimate update
P(:,:,i+1)=P(:,:,i+1)-K*H*P(:,:,i+1);
%Tracking error covariance update
P_updated(:,:,i+1)=P(:,:,i+1);
end
%====================================================================
%%Estimate error
x_error=x-x_hat;
%%Graph1practical and tracking position figure(1)
plot(x(1,:),x(3,:),'r');
hold on;
plot(x_hat(1,:),x_hat(3,:),'g.');
title('2D Target Position')
legend('Practical Position','Tracking Position') xlabel('X axis[m]')
ylabel('Y axis[m]')
hold off;
%%Graph2
figure(2)
plot(t,x(1,:)),grid on;
hold on;
plot(t,x_hat(1,:),'r'),grid on;
title('Practical and Tracking Position on X axis') legend('Practical Position','Tracking Position') xlabel('Time[sec]')
ylabel('Position[m]')
hold off;
%%Graph3
figure(3)
plot(t,x_error(1,:)),grid on;
title('Position Error on X axis')
xlabel('Time[sec]')
ylabel('Position RMSE[m]')
hold off;
%%Graph4
figure(4)
plot(t,x(2,:)),grid on;
hold on;
plot(t,x_hat(2,:),'r'),grid on;
title('Practical and Tracking Velocity on X axis') legend('Practical Velocity','Tracking Velocity') xlabel('Time[sec]')
ylabel('V elocity[m/sec]')
%%Graph5
figure(5)
plot(t,x_error(2,:)),grid on;
title('Velocity Error on X axis')
xlabel('Time[sec]')
ylabel('V elocity RMSE[m/sec]')
hold off;
%==================================================================== ======
一个应用实例详解卡尔曼滤波及其算法实现 标签：算法filtermatlabalgorithm优化工作 2012-05-14 10:48 75511人阅读评论(25) 收藏举报分类： 数据结构及其算法（4） 为了可以更加容易的理解卡尔曼滤波器，这里会应用形象的描述方法来讲解，而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是，他的5条公式是其核心内容。结合现代的计算机，其实卡尔曼的程序相当的简单，只要你理解了他的那5条公式。 在介绍他的5条公式之前，先让我们来根据下面的例子一步一步的探索。 假设我们要研究的对象是一个房间的温度。根据你的经验判断，这个房间的温度是恒定的，也就是下一分钟的温度等于现在这一分钟的温度（假设我们用一分钟来做时间单位）。假设你对你的经验不是100%的相信，可能会有上下偏差几度。 我们把这些偏差看成是高斯白噪声（White Gaussian Noise），也就是这些偏差跟前后时间是没有关系的而且符合高斯分配（Gaussian Distribution）。另外，我们在房间里放一个温度计，但是这个温度计也不准确的，测量值会比实际值偏差。我们也把这些偏差看成是高斯白噪声。 好了，现在对于某一分钟我们有两个有关于该房间的温度值：你根据经验的预测值（系统的预测值）和温度计的值（测量值）。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。 假如我们要估算k时刻的是实际温度值。首先你要根据k-1时刻的温度值，来预测k时刻的温度。因为你相信温度是恒定的，所以你会得到k时刻的温度预测值是跟k-1时刻一样的，假设是23度，同时该值的高斯噪声的偏差是5度（5是这样得到的：如果k-1时刻估算出的最优温度值的偏差是3，你对自己预测的不确定度是4度，他们平方相加再开方，就是5）。然后，你从温度计那里得到了k时刻的温度值，假设是25度，同时该值的偏差是4度。 由于我们用于估算k时刻的实际温度有两个温度值，分别是23 度和25度。究竟实际温度是多少呢？相信自己还是相信温度计呢？究竟相信谁多一点，我们可以用他们的covariance（协方差）来判断。因为Kg^2=5^2/(5^2+4^2)，所以Kg=0.78，我们可以估算出k时刻的实际温度值是：23+0.78*(25-23)=24.56度。 可以看出，因为温度计的covariance比较小（比较相信温度计），所以估算出的最优温度值偏向温度计的值。 现在我们已经得到k时刻的最优温度值了，下一步就是要进入k+1时刻，进行新的最优估算。到现在为止，好像还没看到什么自回归的东西出现。对了，在进入k+1时刻之前，我们还要算出k时刻那个最优值（24.56 度）的偏差。算法如下：((1-Kg)*5^2)^0.5=2.35。这里的5就是上面的k时刻你预测的那个23度
IIR 低通滤波器的设计程序为： Ft=8000; Fp=1000; Fs=1200; As=100 ; Ap=1; wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft; fp=2*Fp*tan(wp/2); fs=2*Fs*tan(ws/2); [n11,wn11]=buttord(wp,ws,1,50,'s'); [b11,a11]=butter(n11,wn11,'s'); [num11,den11]=bilinear(b11,a11,0.5); [h,w]=freqz(num11,den11); axes(handles.axes1); plot(w*8000*0.5/pi,abs(h)); xlabel('Frequency/Hz'); ylabel('Magnitude'); title('巴特沃斯数字低通滤波器'); 巴特沃斯带通滤波器设计程序为: Ft=8000; Fp1=1200; Fp2=3000; Fs1=1000; Fs2=3200; As=100; Ap=1; wp1=tan(pi*Fp1/Ft); wp2=tan(pi*Fp2/Ft); ws1=tan(pi*Fs1/Ft); ws2=tan(pi*Fs2/Ft); w=wp1*wp2/ws2;
bw=wp2-wp1; wp=1; ws=(wp1*wp2-w.^2)/(bw*w); [n12,wn12]=buttord(wp,ws,1,50,'s'); [b12,a12]=butter(n12,wn12,'s'); [num2,den2]=lp2bp(b12,a12,sqrt(wp1*wp2),bw); [num12,den12]=bilinear(num2,den2,0.5); [h,w]=freqz(num12,den12); plot(w*8000*0.5/pi,abs(h)); axis([0 4000 0 1.5]); xlabel('Frequency/Hz'); ylabel('Magnitude'); title('巴特沃斯数字带通滤波器'); IIR 高通滤波器的设计程序为： Ft=8000; Fp=4000; Fs=3500; wp1=tan(pi*Fp/Ft); ws1=tan(pi*Fs/Ft); wp=1; ws=wp1*wp/ws1; [n13,wn13]=cheb1ord(wp,ws,1,50,'s'); [b13,a13]=cheby1(n13,1,wn13,'s'); [num,den]=lp2hp(b13,a13,wn13); [num13,den13]=bilinear(num,den,0.5); [h,w]=freqz(num13,den13); axes(handles.axes1); plot(w*21000*0.5/pi,abs(h)); xlabel('Frequency/Hz'); ylabel('Magnitude'); title('切比雪夫Ⅰ型数字高通滤波器');
中值滤波原理及MATLAB实现 摘要：图像是一种重要的信息源，通过图像处理可以帮助人们了解信息的内涵。本文将纯净的图像加入椒盐噪声,然后采用中值滤波的方法对其进行去噪。中值滤波是一种常用的非线性信号处理技术，在图像处理中，它对滤除脉冲干扰噪声最为有效。文章阐述了中值滤波的原理、算法以及在图像处理中的应用。MATLAB是一种高效的工程计算语言，在数据处理、图像处理、神经网络、小波分析等方面都有广泛的应用。 关键词：图像,中值滤波,去噪,MATLAB 1. 引言 20世纪20年代,图像处理首次得到应用。上个世纪60年代中期,随着计算机科学的发展和计算机的普及,图像处理得到广泛的应用。60年代末期,图像处理技术不断完善,逐渐成为一个新兴的学科。图像处理中输入的是质量低的图像,输出的是改善质量后的图像。 为了改善图像质量,从图像中提取有效信息,必须对图像进行去噪预处理。根据噪声频谱分布的规律和统计特征以及图像的特点,出现了多种多样的去噪方法。经典的去噪方法有：空域合成法,频域合成法和最优合成法等,与之适应的出现了许多应用方法,如均值滤波器,中值滤波器,低通滤波器,维纳滤波器,最小失真法等。这些方法的广泛应用,促进数字信号处理的极大发展,显著提高了图像质量。 2. 中值滤波 在图像滤波中，常用的方法是线性滤波技术和非线性滤波技术，线性滤波以其完美的理论基础，数学处理简单、易于采用和硬件实现等优点，一直在图像滤波领域中占有重要的地位。线性滤波对加性高斯噪声有较好的平滑作用，但对脉冲信号和其它形式的高频分量抑制效果较差，且模糊信号边缘。非线性滤波是基于对输入信号序列的一种非线性投影关系，常把某一特定的噪声近似为零而保留信号的重要特征，一定程度上克服线性滤波器的不足，非线性滤波早期运用较多的是中值滤波器，其应用于多维信号处理时，对窄脉冲信号具有良好的抑制能力，但
《数字信号处理课程设计报告》 课程设计题目：数字带通滤波设计 学院：信息工程学院 专业：通信工程 班级： 学生姓名: 指导老师： 日期：2012年5月4日至17日
目录 1. 课程设计的目的和意义 (3) 1.1课程目的 (3) 1.2实验意义 (3) 2.课程设计题目描述及要求 (4) 3.实验内容 (4) 3.1数字滤波器的简介 (4) 3.2 椭圆数字带通滤波器设计原理 (4) 3.3实验流程框图 (5) 3.4实验步骤 (7) 3.5 程序代码 (8) 3.6实验仿真结果图 (9) 3.7实验结果分析 (10) 4.实验总结 (10) 参考文献 (11)
1.课程设计的目的和意义 1.1课程目的 (1) 学习椭圆模拟带通滤波器的设计。 (2) 使用椭圆模拟带通滤波器逼近数字带通滤波器。 (3) 双线性变换法的应用原理。 (4) 使用matlab对滤波器进行仿真和频谱分析。 1.2实验意义 数字滤波器是个离散时间系统（按预定的算法，将输入离散时间信号换为所要求的输出离散时间信号的特定功能装置）。应用数字滤波器处理模拟信号时，首先须对输入模拟信号进行限带、抽样和模数转换。数字滤波器输入信号的抽样率应大于被处理信号带宽的两倍，其频率响应具有以抽样频率为间隔的周期重复特性，且以折叠频率即1/2抽样频率点呈镜像对称。为得到模拟信号，数字滤波器处理的输出数字信号须经数模转换、平滑。数字滤波器具有高精度、高可靠性、可程控改变特性或复用、便于集成等优点。数字滤波器在语言信号处理、图像信号处理、医学生物信号处理以及其他应用领域都得到了广泛应用。 利用Matlab设计滤波器,可以按照设计要求非常方便地调整设计参数,极大地减轻了设计的工作量,有利于滤波器设计的最优化。Matlab因其强大的数据处理功能被广泛应用于工程计算,其丰富的工具箱为工程计算提供了便利,利用Matlab信号处理工具箱可以快速有效地设计各种数字滤波器,设计简单方便。本文介绍了在Matlab环境下滤波器设计的方法和步骤。
毕业设计(论文) UNDERGRADUATE PROJECT (THESIS) 题目: 冲击测试常用滤波算法研究 学院 专业 学号 学生姓名 指导教师 起讫日期
目录 摘要 (2) ABSTRACT (3) 第一章绪论 (4) 1.1课题背景 (4) 1.2国内外相关领域的研究 (4) 1.3主要研究内容与创新 (5) 1.3.1研究内容与意义 (5) 1.3.2课题的创新点 (5) 1.3.3 研究目的与技术指标 (6) 第二章数字滤波基础 (7) 2.1数字滤波算法概念 (7) 2.2数据采样与频谱分析原理 (8) 2.2.1 时域抽样定理 (8) 2.2.2 离散傅立叶变换(DFT) (8) 2.2.3 快速傅立叶变换(FFT) (9) 2.2.4 频谱分析原理 (9) 2.3常用数字滤波算法基础 (10) 2.3.1常用数字滤波算法分类 (10) 2.3.2常用数字滤波算法特点 (11) 2.3.3常用滤波算法相关原理 (13) 2.4 冲击测试采样数据 (16) 2.4.1噪声的特点与分类 (16) 2.4.2冲击测试采样数据特点 (17) 2.5 MATLAB简介 (17) 2.5.1 MATLAB功能简介 (18) 2.5.2 MATLAB的发展 (18) 第三章、冲击测试滤波算法设计及滤波效果分析 (20) 3.1 冲击测试采样数据的分析 (20) 3.2 滤波算法设计及效果分析 (21) 3.2.1 中位值平均法的设计 (21) 3.2.2限幅法和限速法的设计 (23) 3.2.3一阶滞后法的设计 (25) 3.2.4低通法的设计 (26) 第四章结论与展望 (34) 4.1冲击测试的滤波算法总结 (34) 4.2冲击测试的滤波算法展望 (34) 致谢 (36) 参考文献 (37) 附录：程序代码清单 (38)
文件一 % THIS PROGRAM IS FOR IMPLEMENTATION OF DISCRETE TIME PROCESS EXTENDED KALMAN FILTER % FOR GAUSSIAN AND LINEAR STOCHASTIC DIFFERENCE EQUATION. % By (R.C.R.C.R),SPLABS,MPL. % (17 JULY 2005). % Help by Aarthi Nadarajan is acknowledged. % (drawback of EKF is when nonlinearity is high, we can extend the % approximation taking additional terms in Taylor's series). clc; close all; clear all; Xint_v = [1; 0; 0; 0; 0]; wk = [1 0 0 0 0]; vk = [1 0 0 0 0]; for ii = 1:1:length(Xint_v) Ap(ii) = Xint_v(ii)*2; W(ii) = 0; H(ii) = ‐sin(Xint_v(ii)); V(ii) = 0; Wk(ii) = 0; end Uk = randn(1,200); Qu = cov(Uk); Vk = randn(1,200); Qv = cov(Vk); C = [1 0 0 0 0]; n = 100; [YY XX] = EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V); for it = 1:1:length(XX) MSE(it) = YY(it) ‐ XX(it); end tt = 1:1:length(XX); figure(1); subplot(211); plot(XX); title('ORIGINAL SIGNAL'); subplot(212); plot(YY); title('ESTIMATED SIGNAL'); figure(2); plot(tt,XX,tt,YY); title('Combined plot'); legend('original','estimated'); figure(3); plot(MSE.^2); title('Mean square error'); 子文件：：function [YY,XX] = EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V); Ap(2,:) = 0; for ii = 1:1:length(Ap)‐1 Ap(ii+1,ii) = 1;
实验4 基于MATLAB 的FIR 数字滤波器设计 实验目的：加深对数字滤波器的常用指标和设计过程的理解。 实验原理：低通滤波器的常用指标： P P P for H Ω≤Ω+≤Ω≤-,1)(1δδ πδ≤Ω≤Ω≤ΩS S for H ,)( 通带边缘频率P Ω，阻带边缘频率S Ω ，通带起伏 P δ， 通带峰值起伏] )[1(log 2010dB p p δα--=， 阻带起伏s δ，最小阻带衰减])[(log 2010dB s S δα-=。 数字滤波器有IIR 和FIR 两种类型，它们的特点和设计方法不同。 在MATLAB 中，可以用b=fir1(N,Wn,’ftype’,taper) 等函数辅助设计FIR 数字滤波器。N 代表滤波器阶数；Wn 代表滤波器的截止频率(归一化频率)，当设计带通和带阻滤波器时，Wn 为双元素相量；ftype 代表滤波器类型，如’high ’高通，’stop ’带阻等；taper 为窗函数类型，默认为海明窗，窗系数需要实现用窗函数blackman, hamming,hanning chebwin, kaiser 产生。 例1 用凯塞窗设计一FIR 低通滤波器，通带边界频率π3.0=Ωp ，阻带边界频率π5.0=Ωs ，阻带衰减 不小于50dB 。 解 首先由过渡带宽和阻带衰减 来决定凯塞窗的N 和 π2.0=Ω-Ω=?Ωp s ， ， S P P S Passband Stopband Transition band Fig 1 Typical magnitude specification for a digital LPF
卡尔曼滤波器及其简matlab仿真
卡尔曼滤波器及其简matlab仿真 一、卡尔曼滤波的起源 谈到信号的分析与处理，就离不开滤波两个字。通常，信号的频谱处于有限的频率范围内，而噪声的频谱则散布在很广的频率范围内，为了消除噪声，可以进行频域滤波。但在许多应用场合，需要直接进行时域滤波，从带噪声的信号中提取有用信号。虽然这样的过程其实也算是对信号的滤波，但其所依据的理论，即针对随机信号的估计理论，是自成体系的。人们对于随机信号干扰下的有用信号不能“确知”，只能“估计”。为了“估计”，要事先确定某种准则以评定估计的好坏程度。 1960年卡尔曼发表了用递归方法解决离散数据线性滤波问题的论文A New Approach to Linear Filtering and Prediction Problems （线性滤波与预测问题的新方法），在这篇文章里一种克服了维纳滤波缺点的新方法被提出来，这就是我们今天称之为卡尔曼滤波的方法。卡尔曼滤波应用广泛且功能强大，它可以估计信号的过去和当前状态甚至能估计将来的状态即使并不知道模型的确切性质。 其基本思想是以最小均方误差为最佳估计准则，采用信号与噪声的状态空间模型利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计，求出当前时刻的估计值。算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的估计。 对于解决很大部分的问题，它是最优，效率最高甚至是最有用的。它的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 卡尔曼滤波不要求保存过去的测量数据，当新的数据到来时，根据新的数据和前一时刻的储值的估计，借助于系统本身的状态转移方程，按照一套递推公式，即可算出新的估值。卡尔曼递推算法大大减少了滤波装置的存储量和计算量，并且突破了平稳随机过程的限制，使卡尔曼滤波器适用于对时变信号的实时处理。
matlab实现中值滤波去除脉冲噪声matlab小程序（图像处理）2010-04-1612:58:44阅读8评论0字号：大中小 实验原理：中值滤波器是将领域内像素灰度的中值代替该像素的值,对处理脉冲噪声(椒盐噪声)非常有效。为了对一幅图像上的某个点进行中值滤波处理，必须先将掩模内欲求的像素及其领域的像素值排序，确定出中值，主要功能是使拥有不同灰度的点看起来更接近于它的邻近值。 程序说明：函数名为mid(pic_name,s)的函数，其中参数pic_name为读入的图像，s为掩模矩阵的边长，由用户自行决定。 实验说明：随着掩模矩阵的变大，我们可以看到脉冲噪声去除得更加理想，但同时图像会变得更模糊，因为各点像素与其邻域更为接近，因此，进行中值滤波时选择一个适合的掩模矩阵十分重要。另外，我们看到图像的边界处出现了黑色的斑点，这是由于我采用了0来直译边界，这种影响可用镜像反射方式对称地沿其边界扩展来减弱。 另附：其实本实验可以完全由matlab中的函数median或medfilt2简单实现，此处写出内部处理过程，主要是为了让大家理解中值滤波的具体处理过程。 程序源代码： function mid(pic_name,s) close all; s=double(s); X=imread(pic_name); Y1=imnoise(X,'salt&pepper',0.2);%对读入的图像加脉冲噪声 figure; imshow(uint8(Y1)); Y1=double(Y1); [m,n]=size(X); s2=round(s/2); s3=round(s*s/2);%中值像素点的位置
卡尔曼滤波器及其简matlab仿真 一、卡尔曼滤波的起源 谈到信号的分析与处理，就离不开滤波两个字。通常，信号的频谱处于有限的频率范围内，而噪声的频谱则散布在很广的频率范围内，为了消除噪声，可以进行频域滤波。但在许多应用场合，需要直接进行时域滤波，从带噪声的信号中提取有用信号。虽然这样的过程其实也算是对信号的滤波，但其所依据的理论，即针对随机信号的估计理论，是自成体系的。人们对于随机信号干扰下的有用信号不能“确知”，只能“估计”。为了“估计”，要事先确定某种准则以评定估计的好坏程度。 1960年卡尔曼发表了用递归方法解决离散数据线性滤波问题的论文A New Approach to Linear Filtering and Prediction Problems（线性滤波与预测问题的新方法），在这篇文章里一种克服了维纳滤波缺点的新方法被提出来，这就是我们今天称之为卡尔曼滤波的方法。卡尔曼滤波应用广泛且功能强大，它可以估计信号的过去和当前状态甚至能估计将来的状态即使并不知道模型的确切性质。 其基本思想是以最小均方误差为最佳估计准则，采用信号与噪声的状态空间模型利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计，求出当前时刻的估计值。算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的估计。 对于解决很大部分的问题，它是最优，效率最高甚至是最有用的。它的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 卡尔曼滤波不要求保存过去的测量数据，当新的数据到来时，根据新的数据和前一时刻的储值的估计，借助于系统本身的状态转移方程，按照一套递推公式，即可算出新的估值。卡尔曼递推算法大大减少了滤波装置的存储量和计算量，并且突破了平稳随机过程的限制，使卡尔曼滤波器适用于对时变信号的实时处理。 二、卡尔曼滤波的原理
clear all v=150; %%目标速度 v_sensor=0;%%传感器速度 t=1; %%扫描周期 xradarpositon=0; %%传感器坐标yradarpositon=0; %% ppred=zeros(4,4); Pzz=zeros(2,2); Pxx=zeros(4,2); xpred=zeros(4,1); ypred=zeros(2,1); sumx=0; sumy=0; sumxukf=0; sumyukf=0; sumxekf=0; sumyekf=0; %%%统计的初值 L=4; alpha=1; kalpha=0; belta=2; ramda=3-L; azimutherror=0.015; %%方位均方误差rangeerror=100; %%距离均方误差processnoise=1; %%过程噪声均方差 tao=[t^3/3 t^2/2 0 0; t^2/2 t 0 0; 0 0 t^3/3 t^2/2; 0 0 t^2/2 t]; %% the input matrix of process G=[t^2/2 0 t 0 0 t^2/2 0 t ]; a=35*pi/180; a_v=5/100; a_sensor=45*pi/180; x(1)=8000; %%初始位置
y(1)=12000; for i=1:200 x(i+1)=x(i)+v*cos(a)*t; y(i+1)=y(i)+v*sin(a)*t; end for i=1:200 xradarpositon=0; yradarpositon=0; Zmeasure(1,i)=atan((y(i)-yradarpositon)/(x(i)-xradarpositon))+random('Normal',0,azimutherror,1,1); Zmeasure(2,i)=sqrt((y(i)-yradarpositon)^2+(x(i)-xradarpositon)^2)+random('Normal',0,rangeerror,1,1); xx(i)=Zmeasure(2,i)*cos(Zmeasure(1,i));%%观测值 yy(i)=Zmeasure(2,i)*sin(Zmeasure(1,i)); measureerror=[azimutherror^2 0;0 rangeerror^2]; processerror=tao*processnoise; vNoise = size(processerror,1); wNoise = size(measureerror,1); A=[1 t 0 0; 0 1 0 0; 0 0 1 t; 0 0 0 1]; Anoise=size(A,1); for j=1:2*L+1 Wm(j)=1/(2*(L+ramda)); Wc(j)=1/(2*(L+ramda)); end Wm(1)=ramda/(L+ramda); Wc(1)=ramda/(L+ramda);%+1-alpha^2+belta; %%%权值 if i==1 xerror=rangeerror^2*cos(Zmeasure(1,i))^2+Zmeasure(2,i)^2*azimutherror^2*sin(Zmeasure(1,i))^2; yerror=rangeerror^2*sin(Zmeasure(1,i))^2+Zmeasure(2,i)^2*azimutherror^2*cos(Zmeasure(1,i))^2; xyerror=(rangeerror^2-Zmeasure(2,i)^2*azimutherror^2)*sin(Zmeasure(1,i))*cos(Zmeasure(1,i)); P=[xerror xerror/t xyerror xyerror/t; xerror/t 2*xerror/(t^2) xyerror/t 2*xyerror/(t^2); xyerror xyerror/t yerror yerror/t;
数字图像处理实验二 15生医 一、实验内容 产生教材104页题图4.18（右图）所示的二值图像（白为1，黑为0），编程实现习题4.18所要求的处理（3x3的平均滤波和中值滤波）功能（图像四周边界不考虑，处理结果按四舍五入仍取0或1），显示处理前后的图像，比较其异同。 二、运行环境 MATLAB R2014a 三、运行结果及分析 1. 四种不同的窗的3x3平均滤波 ①在 MATLAB 图形窗界面进行放大可以看出四者之间的差别： 4领域与8邻域之间没有明显差别，但是加权与未加权之间的差别较为明显，体现在： 加权后每个矩形块的四个尖角部分都被保留了下来4邻域平均滤波后图 像8邻域平均滤波后图像 4邻域加权平均滤波后图像8邻域加权平均滤波后图像
（图像四周边界不考虑），而未加权的尖角处黑色变为白色。 ②原因分析： 加权后尖角处原来白色的点（1）进行计算3/5=0.6四舍五入后值为1，保持白色，原来黑色的点（0）进行计算2/5=0.4四舍五入后值为0，保持黑色；而未加权尖角处无论原来是黑色还是白色，进行计算 2/4=0.5四舍五入后值为1，所以原先的黑色（0）也变成了白色（1）。 ③下图为放大后的截图： 2.中值滤波与原图像的对比
①在 MATLAB图形窗界面进行放大后可观察出： 使用3x3 方形中值滤波模板的效果与4领域、8领域加权平均滤波的 效果相同，每个矩形块的四个尖角部分都被保留了下来（图像四周边界不考虑）。 ②原因分析： 套用3x3方形中值滤波模板后，尖角处原来白色的点（1）在窗内1多于0，取中值后仍保持白色，原来黑色的点（0）在窗内0多于1，取中值后仍保持白色。 ③下图为放大后的截图： 原图像中值滤波后图像
基于matlab的fir数字滤波器设计 sea （湖南涉外经济学院湖南长沙 421000） 摘要：介绍了应用Matlab语言设计FIR数字滤波器时采用直接程序设计法。同时介绍了FIR数字滤波器几种设计方法的函数调用格式；通过实例，给出了程序设计法的详细步骤，并在Matlab的Simulink环境下，对所设计的滤波器进行了仿真。 关键词：数字滤波器；Matlab Abstract: It Introduced the application of Matlab language when designing FIR digital filter with direct procedures. At the same time it introduced several formats of function calls when designing FIR digital filter; through some examples of programming method it gives the detail steps of Simulink in Matlab environment, and to simulate the designed filter. Keywords: digital filter MATLAB 引言 1. MATLAB的特点 MATLAB(Matrix Laboratory)机矩阵实验室是由美国MathWorks公司推出的一款集数值分析、矩阵计算、信号处理和图形显示于一体的工具软件。它的命令语句功能非常强大,包含了大量高度集成的可直接调用的函数，高效简洁；另一方面，它又是一个开放系统，针对不同的学科，推出了不同的工具箱。自1984年推向市场以来，经过十几年的发展和竞争，现已成为国际认可(IEEE)的最优化的科技应用软件一。 正是由于MATLAB具有良好的扩展性以及强大的数据分析和处理能力，现已广泛应用于矩阵代数、数值计算、数字信号处理、振动理论、神经网络控制、动态仿真等领域。MATLAB是用M语言编程，尽管不能在M文件中直接调用C语言程序，但可以通过MATLAB提供的应用编程接口(API)来与外部接口，在MATLAB环境中调用C语言或Fortran程序、输入或输出数据以及与其他软件程序建立客户/服务 器关系。在MATLAB中调用C语言程序，必须通过MEX文件来实现。
function g=QuiMedFil2(f) g=f; t=2; [m,n]=size(f); for i=2:m-1 for j=2:n-1 if(j==2) A=f(i-1:i+1,j-1:j+1); hist=imhist(A); mdn=median(A(:)); %mdn=Med(f(i-1:i+1,j-1:j+1))+1; L=find(f(i-1:i+1,j-1:j+1)
while (Ltmdn+hist(mdn+1))<=th Ltmdn=Ltmdn+hist(mdn+1); mdn=mdn+1; end g(i,j)=mdn; end end t=t+1;%t-line; end ======================================= =========================实现过程如下 rgb=imread('Figure1.JPG'); %此处通过matlab读入任何一幅图像 >> I=rgb2gray(rgb); %将图像灰度化 >> I=imresize(I,0.1); %若图像过大可进行适当调整，或省略此步 >> g=QuiMedFil2(I) %调用所编程序执行中值滤波快速算法 figure，imshow(I),figure,imshow(g) %显示原图像和滤波后的图像 =============================================================================== =================================================== 与传统中值滤波算法的同窗口的计算时间对比 快速算法略高于一般算法 =============================================================================== =====================================================菱形5*5中值滤波程序function g=MedFilRho(f) g=f; g=[]; [m,n]=size(f); for i=3:(m-2) for j=3:(n-2) B=[f(i,j),f(i-2,j),f(i-1,j),f(i+1,j),f(i+2,j),f(i,j-1),f(i,j-2),f(i,j+1),f(i,j+2),f(i-1,j-1),f(i+1,j+1),f(i-1,j+1),f(i+1,j-1)]; g(i,j)=median(B); end end 备注：程序为原创，原文地址： https://www.doczj.com/doc/9f12250034.html,/walqxlqxaw/blog/item/aabb517f877e6f190dd7daad.ht ml
目录 1 引言 (1) 2 设计任务 (2) 2.1设计内容 (2) 2.2设计要求 (2) 3 语音信号的采集及时频分析 (3) 3.1语音信号的采集 (3) 3.2语音信号的时频分析 (3) 4 基于MATLAB的数字滤波器的设计 (5) 4.1数字滤波器的设计 (5) 4.1.1数字滤波器的基本概念 (5) 4.1.2 IIR滤波器设计思想 (5) 4.2IIR数字滤波器设计 (5) 4.2.1 IIR低通滤波器设计 (5) 4.2.2 IIR带通滤波器设计 (7) 4.2.3 IIR带通滤波器设计 (9) 5 合成信号及其滤波 (12) 5.1合成信号 (12) 5.2合成信号滤波 (13) 6 设计系统界面 (15) 6.1系统界面设计工具—GUI概述 (15) 6.2界面设计及使用说明 (15) 7 心得体会 (18) 参考文献 (19) 附录 (20)
1 引言 数字信号处理是利用计算机或专用处理设备，以数值计算的方法对信号进行采集、抽样、变换、综合、估值与识别等加工处理，借以达到提取信息和便于应用的目的。它在语音、雷达、图像、系统控制、通信、航空航天、生物医学等众多领域都获得了极其广泛的应用。具有灵活、精确、抗干扰强、设备尺寸小、造价低、速度快等优点。 数字滤波器, 是数字信号处理中及其重要的一部分。随着信息时代和数字技术的发展，受到人们越来越多的重视。数字滤波器可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。 IIR滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。 MATALB 可以创建图形用户界面GUI (GraphicalUser Interface) ,它是用户和计算机之间交流的工具。MATLAB 将所有GUl 支持的用户控件都集成在这个环境中并提供界面外观、属性和行为响应方式的设置方法,随着版本的提高,这种能力还会不断加强。而且具有强大的绘图功能,可以轻松的获得更高质量的曲线图。 滤波器的设计可以通过软件或设计专用的硬件两种方式来实现。随着MATLAB软件及信号处理工具箱的不断完善，MATLAB很快成为应用学科等领域不可或缺的基础软件。它可以快速有效地实现数字滤波器的设计、分析和仿真，极大地减轻了工作量,有利于滤波器设计的最优化。
MATLAB 对卡尔曼滤波器的仿真实现 刘丹，朱毅，刘冰 武汉理工大学信息工程学院，武汉（430070） E-mail ：liudan_ina@https://www.doczj.com/doc/9f12250034.html, 摘 要：本文以卡尔曼滤波器原理为理论基础，用MATLAB 进行卡尔曼滤波器仿真、对比卡尔曼滤波器的预测效果，对影响滤波其效果的各方面原因进行讨论和比较，按照理论模型进行仿真编程，清晰地表述了编程过程。 关键词：数字信号处理；卡尔曼滤波器；MATLAB ；仿真过程 中图分类号: TN912.3 1. 引言 随着信息时代和数字世界的到来，数字信号处理已成为当今一门极其重要的学科和技术领域。数字信号处理已在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理中,数字滤波占有极其重要的地位,目前对数字滤波器的设计有多种方法,其中著名的MATLAB 软件包在多个研究领域都有着广泛的应用,它的频谱分析[1]和滤波器的分析设计功能很强,从而使数字信号处理变得十分简单、直观。本文分析了数字滤波器的设计方法,举出了基于MATLAB 软件的信号处理工具在数字滤波器设计中的应用。 2. 卡尔曼滤波基本原理 卡尔曼滤波过程实际上是获取维纳解的递推运算过程[2]。从维纳解导出的卡尔曼滤波器实际上是卡尔曼滤波过程结束后达到稳态的情况，这时Kalman Filtering 的结果与Wiener Solution 是相同的[3]。具体推导如下： )()1|1(?)|(?n Gy n n x f n n x +??= )|(?)()(n n x n x n e ?= 已知由此求c a cG a f F G n e E n ,)1(( ..min )]([)(2?=??→?==ε 由 f G f G ,0??????????=??εε ⑴ )]1|1(?)()[()1|1(?)|(????+??=n n x ac n y n G n n x a n n x 可以是时变的，非平稳的随机信号 ⑵ Q n a n P +?=)1()(2 ε均为正数。 ⑶ ) () ()(2n P C R n CP n G += ⑷ )()](1[)()(n P n CG n G C P n ??== ε )(n G 是个随时间变化的量，每次输入输出，)(n G 就调整一次，并逐渐逼近Kalman Filter 的增益G ，而)1()(? 实验7\8基于MATLAB勺数字滤波器设计实验目的：加深对数字滤波器的常用指标和设计过程的理解。 实验原理：低通滤波器的常用指标： 1 一6P 兰G(e^) ≤ 1 + 6P , for 国≤ ωP G(J") ≤ 6s, for 国s ≤ ⑷≤ ∏ 通带边缘频率：'P ,阻带边缘频率：'s， 通带起伏：J P，通带峰值起伏： C(P= —20 IOg io (^-OP )【d B 】阻带起伏.冠S PaSSband StOPband Tran Siti on band Fig 7.1 TyPiCaI magn itude SPeCifiCati On for a digital LPF :S = -20 log ιo(r)[dB 】 O 数字滤波器有IIR和FlR两种类型，它们的特点和设计方法不同。 在MATLAB^,可以用［b , a］=butter ( N,Wr)等函数辅助设计IIR数字滤波器，也可以用b=fir1(N,Wn, 'type ')等函数辅助设计FIR数字滤波器。 实验内容：利用MATLAB编程设计一个数字带通滤波器，指标要求如下： 通带边缘频率：??P1=0.45^,?? P2=0?65 二，通带峰值起伏：［dB】O 阻带边缘频率：'s1 0.3…,'s2 0.75…，最小阻带衰减：-S 4°[dB] O 分别用IIR和FlR两种数字滤波器类型进行设计。 实验要求：给出IIR数字滤波器参数和FIR数字滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。 实验内容： IRR代码： wp=[0.45*pi,0.65*pi]; ws=[0.3*pi,0.75*pi]; Ap=1; A S=40; [N,Wc]=buttord(wp∕pi,ws∕pi,Ap,As); [b,a]=butter(N,Wc)%[b,a] = butter( n, Wn,'ftype') 最小阻带衰减: 1线性平滑滤波器 用MA TLAB实现领域平均法抑制噪声程序： I=imread(' c4.jpg '); subplot(231) imshow(I) title('原始图像') I=rgb2gray(I); I1=imnoise(I,'salt & pepper',0.02); subplot(232) imshow(I1) title(' 添加椒盐噪声的图像') k1=filter2(fspecial('average',3),I1)/255; %进行3*3模板平滑滤波 k2=filter2(fspecial('average',5),I1)/255; %进行5*5模板平滑滤波 k3=filter2(fspecial('average',7),I1)/255; %进行7*7模板平滑滤波 k4=filter2(fspecial('average',9),I1)/255; %进行9*9模板平滑滤波 subplot(233),imshow(k1);title('3*3 模板平滑滤波'); subplot(234),imshow(k2);title('5*5 模板平滑滤波'); subplot(235),imshow(k3);title('7*7 模板平滑滤波'); subplot(236),imshow(k4);title('9*9 模板平滑滤波'); 2.中值滤波器 用MA TLAB实现中值滤波程序如下： I=imread(' c4.jpg '); I=rgb2gray(I); J=imnoise(I,'salt&pepper',0.02); subplot(231),imshow(I);title('原图像'); subplot(232),imshow(J);title('添加椒盐噪声图像'); k1=medfilt2(J); %进行3*3模板中值滤波 k2=medfilt2(J,[5,5]); %进行5*5模板中值滤波 k3=medfilt2(J,[7,7]); %进行7*7模板中值滤波 k4=medfilt2(J,[9,9]); %进行9*9模板中值滤波 subplot(233),imshow(k1);title('3*3模板中值滤波'); subplot(234),imshow(k2);title('5*5模板中值滤波'); subplot(235),imshow(k3);title('7*7模板中值滤波'); subplot(236),imshow(k4);title('9*9 模板中值滤波'); 3状态统计滤波器：ordfilt2函数 Y=ordfilt2(X,order,domain) 由domain中非0元素指定邻域的排序集中的第order个元素代替X中的每个元素。Domain 是一个仅包括0和1的矩阵，1仅定义滤波运算的邻域。 Y=ordfilt2(X,order,domain,S) S与domain一样大，用与domain的非0值相应的S的值作为附加补偿。 4二维自适应除噪滤波器：wiener2函数 wiener2函数估计每个像素的局部均值与方差，该函数用法如下：实验报告基于MATLAB的数字滤波器设计
MATLAB滤波程序
相关主题
文本预览