维纳滤波器-matlab实现 (1)

实验报告册

数字图形图像处理

维纳滤波器matlab实现

学院：人民武装学院学院

专业：计算机科学与技术

班级： 11级计科班

学号： 1120070544 学生姓名：苏靖

指导教师：

维纳滤波的原理及其matlab 实现，以案例的形式展示FIR 维纳滤波的特性。

2.维纳滤波概述

维纳（Wiener ）是用来解决从噪声中提取信号的一种过滤（或滤波）方法。这种线性滤波问题，可以看做是一种估计问题或一种线性估计问题。

一个线性系统，如果它的单位样本响应为)(n h ，当输入一个随机信号)(n x ，且

)()()(n v n s n x += （1） 其中)(n x 表示信号，)(n v )表示噪声，则输出)(n y 为

∑-=m

m n x m h n y )()()( （2）

我们希望)(n x 通过线性系统)(n h 后得到的)(n y 尽量接近于)(n s ，因此称)(n y 为)(n s 的估计值，用^

)(n s 表示，即

^)()(n s n y = （3） 则维纳滤波器的输入—输出关系可用下面图1表示。

图1

实际上，式（2）所示的卷积形式可以理解为从当前和过去的观察值)(n x ，)1(-n x ，)2(-n x …)(m n x -，…来估计信号的当前值^)(n s 。因此，用)(n h 进行过滤问题实际上是一种统计估计问题。

一般地，从当前的和过去的观察值)(n x ，)1(-n x ，)2(-n x …估计当前的信号值^

)()(n s n y =成为过滤或滤波；从过去的观察值，估计当前的或者将来的信号值)0)(()(^≥+=N N n s n y 称为外推或预测；从过去的观察值，估计过去的信号值)1)(()(^>-=N N n s n y 称为平滑或内插。因此维纳滤波器又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓的最佳与最优是以最小均方误差为准则的。

如果我们分别以)(n s 与^)(n s 表示信号的真实值与估计值，而用)(n e 表示他们之间的误差，即

)()()(^n s n s n e -= （4） 显然)(n e 可能是正值，也可能是负值，并且它是一个随机变量。因此，用它的均方误差来

表达误差是合理的，所谓均方误差最小即它的平方的统计期望最小：

min )]([)(2

==n E n e ξ （5） 采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。

4.FIR 维纳滤波器的matlab 实现

4.1问题描述

假设一个点目标在x ，y 平面上绕单位圆做圆周运动，由于外界干扰，其运动轨迹发生了偏移。其中，x 方向的干扰为均值为0，方差为0.05的高斯噪声；y 方向干扰为均值为0，方差为0.06的高斯噪声。

1) 产生满足要求的x 方向和y 方向随机噪声500个样本；

2) 明确期望信号和观测信号；

3) 试设计一FIR 维纳滤波器，确定最佳传递函数：1

opt xx xs h R R -=，并用该滤波器处理

观测信号，得到其最佳估计。（注：自行设定误差判定阈值，根据阈值确定滤波器的阶数或传递函数的长度）。

4) 分别绘制出x 方向和y 方向的期望信号、噪声信号、观测信号、滤波后信号、最小均

方误差信号的曲线图；

5) 在同一幅图中绘制出期望信号、观测信号和滤波后点目标的运动轨迹。 4.2 Matlab 仿真及运行结果

用Matlab 实现FIR 滤波器，并将先前随机产生的500个样本输入，得到最佳估计。具体程序如下：

clear;

clf;

sita=0:pi/249.5:2*pi;

xnoise=sqrt(0.05)*randn(1,500);%产生x 轴方向噪声

ynoise=sqrt(0.06)*randn(1,500);%产生y 轴方向噪声

x=cos(sita)+xnoise;%产生x 轴方向观测信号

y=sin(sita)+ynoise;%产生y 轴方向观测信号

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

rxx=xcorr(x);

for i=1:100

for j=1:100

mrxx(i,j)=rxx(500-i+j);

end

end

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

rxd=xcorr(x,xd);

for i=1:100

mrxd(i)=rxd(499+i);

end

hoptx=inv(mrxx)*mrxd';%由维纳-霍夫方程得到的x方向上的滤波器最优解fx=conv(x,hoptx);%滤波后x方向上的输出

nx=sum(abs(xd).^2);

eminx=nx-mrxd*hoptx;%x方向上最小均方误差

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

ryy=xcorr(y);

for i=1:100

for j=1:100

mryy(i,j)=ryy(500-i+j);

end

end

yd=sin(sita);

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

ryd=xcorr(y,yd);

for i=1:100

mryd(i)=ryd(499+i);

end

hopty=inv(mryy)*mryd';%由维纳-霍夫方程得到的y方向上的滤波器最优解fy=conv(y,hopty);%滤波后y方向上的输出

ny=sum(abs(yd).^2);

eminy=ny-mryd*hopty;%y方向上最小均方误差

subplot(2,4,1)

plot(xd);

title('x方向期望信号');

plot(xnoise);

title('x方向噪声信号'); subplot(2,4,3)

plot(x);

title('x方向观测信号'); subplot(2,4,4)

n=0:500;

plot(n,eminx);

title('x方向最小均方误差'); subplot(2,4,5)

plot(yd);

title('y方向期望信号'); subplot(2,4,6)

plot(ynoise);

title('y方向噪声信号'); subplot(2,4,7)

plot(y);

title('y方向观测信号'); subplot(2,4,8)

plot(n,eminy);

title('y方向最小均方误差'); figure;

plot(xd,yd,'k');

hold on;

plot(x,y,'b:');

hold on;

plot(fx,fy,'g-');

title('最终结果');

运行结果如下：

图2

x方向及y方向的期望信号、噪声信号、观测信号以及滤波后的最小均方误差如上图2所示。

图3

滤波后的到的信号与原始信号和噪声信号的对比如上图3所示，滤波后的结果与期望信号还是很接近的，整体上达到了最优滤波的效果。

维纳滤波器的局限

维纳滤波复原法存在着几个实质性的局限。第一，最有标准是基于最小均方误差的且对所有误差等权处理，这个标准在数学上可以接受，但却是个不适合人眼的方式，原因在于人类对复原错误的感知在具有一致灰度和亮度的区域中更为严重，而对于出现在暗的和高梯度区域的误差敏感性差得多。第二，空间可变的退化不能用维纳滤波复原法复原，而这样的退化是常见的。第三，维纳滤波不能处理非平稳信号和噪声。

四、模拟仿真

运行结果

运行程序代码

clear;

I=imread('img_orignal.tif');

figure;

subplot(2,2,1);imshow(I);title('原图像');

[m,n]=size(I);F=fftshift(fft2(I));

k=0.005;

卡尔曼滤波算法与matlab实现

一个应用实例详解卡尔曼滤波及其算法实现 标签：算法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度

几个分形的matlab实现

几个分形得matlab实现 摘要:给出几个分形得实例,并用matlab编程实现方便更好得理解分形,欣赏其带来得数学美感 关键字:Koｃh曲线实验图像 一、问题描述: 从一条直线段开始,将线段中间得三分之一部分用一个等边三角形得两边代替,形成山丘形图形如下 ?图1 在新得图形中,又将图中每一直线段中间得三分之一部分都用一个等边三角形得两条边代替,再次形成新得图形如此迭代,形成Kocｈ分形曲线。 二、算法分析: 考虑由直线段(2个点)产生第一个图形(５个点)得过程。图1中,设与分别为原始直线段得两个端点,现需要在直线段得中间依次插入三个点,,。显然位于线段三分之一处,位于线段三分 之二处,点得位置可瞧成就是由点以点为轴心,逆时针旋转600而得。旋转由正交矩阵 实现。 算法根据初始数据(与点得坐标),产生图1中5个结点得坐标、结点得坐标数组形成一个矩阵,矩阵得第一行为得坐标,第二行为得坐标……,第五行为得坐标。矩阵得第一列元素分别为5个结点得坐标,第二列元素分别为5个结点得坐标。 进一步考虑Kｏｃh曲线形成过程中结点数目得变化规律。设第次迭代产生得结点数为,第次迭代产生得结点数为,则与中间得递推关系为。 三、实验程序及注释: p=［0 0;10 0］; %P为初始两个点得坐标,第一列为ｘ坐标,第二列为y坐标 n=2; %n为结点数 A=[cos(pi/３) —ｓin(pi/3);siｎ(ｐi/３) cos(pi/３)]; ％旋转矩阵 ｆor k＝1:4 ｄ=diｆf(p)/3; ％ｄｉff计算相邻两个点得坐标之差,得到相邻两点确定得向量 %则d就计算出每个向量长度得三分之一,与题中将线段三等分对应 m=4＊n-3; ％迭代公式 q=p(1:ｎ—1,:); %以原点为起点,前n—１个点得坐标为终点形成向量 p(5:4:m,:)=ｐ(2:n,:); %迭代后处于４k+１位置上得点得坐标为迭代前得相应坐标 p(2:4:m,:)=q+d; %用向量方法计算迭代后处于４k+２位置上得点得坐标 p(3:4:m,:)=ｑ+d+ｄ＊A'; %用向量方法计算迭代后处于4ｋ＋3位置上得点得坐标 p(4:４:m,:)＝ｑ＋2*d; %用向量方法计算迭代后处于4k位置上得点得坐标 n=m; %迭代后新得结点数目 end pｌot(ｐ(:,1),p(:,２)) ％绘出每相邻两个点得连线 aｘis([0 10 ０ 10]) 四、实验数据记录: 由第三部分得程序,可得到如下得Ｋoch分形曲线:

简单低通滤波器设计及matlab仿真

东北大学 研究生考试试卷 考试科目： 课程编号： 阅卷人: 考试日期： 姓名：xl 学号： 注意事项 1.考前研究生将上述项目填写清楚. 2.字迹要清楚,保持卷面清洁. 3.交卷时请将本试卷和题签一起上交. 4.课程考试后二周内授课教师完成评卷工作,公共课成绩单与试卷交研究生院培养办公室, 专业课成绩单与试卷交各学院,各学院把成绩单交研究生院培养办公室. 东北大学研究生院培养办公室

数字滤波器设计 技术指标： 通带最大衰减： =3dB ， 通带边界频率： ＝100Hz 阻带最小衰减： ＝20dB 阻带边界频率： ＝200Hz 采样频率：Fs=200Hz 目标： 1、根据性能指标设计一个巴特沃斯低通模拟滤波器。 2、通过双线性变换将该模拟滤波器转变为数字滤波器。 原理： 一、模拟滤波器设计 每一个滤波器的频率范围将直接取决于应用目的，因此必然是千差万别。为了使设计规范化，需要将滤波器的频率参数作归一化处理。设所给的实际频 率为Ω（或f ），归一化后的频率为λ，对低通模拟滤波器令λ＝p ΩΩ/，则1 =p λ， p s s ΩΩ=/λ。令归一化复数变量为p ，λj p =，则p p s j j p Ω=ΩΩ==//λ。所以巴 特沃思模拟低通滤波器的设计可按以下三个步骤来进行。 (1)将实际频率Ω规一化 (2)求Ωc 和N 11010/2-=P C α s p s N λααlg 1 10 110lg 10 /10/--= 这样Ωc 和N 可求。 p x fp s x s f

根据滤波器设计要求＝3dB ，则C ＝1，这样巴特沃思滤波器的设计就只剩一个参数N ，这时 N p N j G 222 )/(11 11)(ΩΩ+= += λλ (3)确定)(s G 因为λj p =，根据上面公式有 N N N p j p p G p G 22)1(11 )/(11)()(-+= += - 由 0)1(12=-+N N p 解得 )221 2exp(πN N k j p k -+=，k ＝1，2， (2) 这样可得 1 )21 2cos(21 ) )((1 )(21+-+-= --= -+πN N k p p p p p p p G k N k k 求得)(p G 后，用p s Ω/代替变量p ，即得实际需要得)(s G 。 二、双线性变换法 双线性变换法是将s 平面压缩变换到某一中介1s 平面的一条横带里，再通过标准变换关系)*1exp(T s z =将此带变换到整个z 平面上去，这样就使s 平面与z 平面之间建立一一对应的单值关系，消除了多值变换性。 为了将s 平面的Ωj 轴压缩到1s 平面的1Ωj 轴上的pi -到pi 一段上，可以通过以下的正切变换来实现： )21 tan(21T T Ω= Ω 这样当1Ω由T pi -经0变化到T pi 时，Ω由∞-经过0变化到∞+，也映射到了整个Ωj 轴。将这个关系延拓到整个s 平面和1s 平面，则可以得到

Matlab实验报告：分形迭代

数学实验报告：分形迭代 练习1 1.实验目的：绘制分形图案并分析其特点。 2.实验内容：绘制Koch曲线、Sierpinski三角形和树木花草图形，观察这些图形的局部和原来分形图形的关系。 3.实验思路：利用函数反复调用自己来模拟分形构造时的迭代过程，当迭代指标n为0时运行作图操作，否则继续迭代。 4.实验步骤： （1）Koch曲线 function koch(p,q,n) % p、q分别为koch曲线的始末复坐标，n为迭代次数 if (n==0) plot([real(p);real(q)],[imag(p);imag(q)]); hold on; axis equal else a=(2*p+q)/3; % 求出从p 到q 的1/3 处端点a b=(p+2*q)/3; % 求出从p 到q 的2/3 处端点b c=a+(b-a)*exp(pi*i/3);% koch(p, a, n-1); % 对pa 线段做下一回合 koch(a, c, n-1); % 对ac 线段做下一回合 koch(c, b, n-1); % 对cb 线段做下一回合 koch(b, q, n-1); % 对bq 线段做下一回合 end （2）Sierpinski三角形 function sierpinski(a,b,c,n) % a、b、c为三角形顶点，n为迭代次数 if (n==0) fill([real(a) real(b) real(c)],[imag(a) imag(b) imag(c)],'b');% 填充三角形abc hold on; axis equal else a1=(b+c)/2; b1=(a+c)/2; c1=(a+b)/2; sierpinski(a,b1,c1,n-1); sierpinski(a1,b,c1,n-1); sierpinski(a1,b1,c,n-1); end （3）树木花草 function grasstree(p,q,n) % p、q分别为树木花草始末复坐标，n为迭代次数

各类滤波器的MATLAB程序清单

各类滤波器的MATLAB程序 一、理想低通滤波器 IA=imread(''); [f1,f2]=freqspace(size(IA),'meshgrid'); Hd=ones(size(IA)); r=sqrt(f1.^2+f2.^2); Hd(r>=0; Y=fft2(double(IA)); Y=fftshift(Y); Ya=Y.*Hd; Ya=ifftshift(Ya); Ia=ifft2(Ya); figure subplot(2,2,1),imshow(uint8(IA)); subplot(2,2,2),imshow(uint8(Ia)); figure surf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 二、理想高通滤波器 IA=imread(''); [f1,f2]=freqspace(size(IA),'meshgrid'); Hd=ones(size(IA)); r=sqrt(f1.^2+f2.^2); Hd(r<=0; Y=fft2(double(IA));

Y=fftshift(Y); Ya=Y.*Hd; Ya=ifftshift(Ya); Ia=real(ifft2(Ya)); figure subplot(2,2,1),imshow(uint8(IA)); subplot(2,2,2),imshow(uint8(Ia)); figure surf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong'); 三、B utterworth低通滤波器 IA=imread(''); [f1,f2]=freqspace(size(IA),'meshgrid'); D=; r=f1.^2+f2.^2; n=4; for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(t^n+1); end end Y=fft2(double(IA)); Y=fftshift(Y); Ya=Y.*Hd; Ya=ifftshift(Ya); Ia=real(ifft2(Ya));

扩展卡尔曼滤波matlab程序

文件一 % 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;

卡尔曼滤波器及其简matlab仿真

卡尔曼滤波器及其简matlab仿真

卡尔曼滤波器及其简matlab仿真 一、卡尔曼滤波的起源 谈到信号的分析与处理，就离不开滤波两个字。通常，信号的频谱处于有限的频率范围内，而噪声的频谱则散布在很广的频率范围内，为了消除噪声，可以进行频域滤波。但在许多应用场合，需要直接进行时域滤波，从带噪声的信号中提取有用信号。虽然这样的过程其实也算是对信号的滤波，但其所依据的理论，即针对随机信号的估计理论，是自成体系的。人们对于随机信号干扰下的有用信号不能“确知”，只能“估计”。为了“估计”，要事先确定某种准则以评定估计的好坏程度。 1960年卡尔曼发表了用递归方法解决离散数据线性滤波问题的论文A New Approach to Linear Filtering and Prediction Problems （线性滤波与预测问题的新方法），在这篇文章里一种克服了维纳滤波缺点的新方法被提出来，这就是我们今天称之为卡尔曼滤波的方法。卡尔曼滤波应用广泛且功能强大，它可以估计信号的过去和当前状态甚至能估计将来的状态即使并不知道模型的确切性质。 其基本思想是以最小均方误差为最佳估计准则，采用信号与噪声的状态空间模型利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计，求出当前时刻的估计值。算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的估计。 对于解决很大部分的问题，它是最优，效率最高甚至是最有用的。它的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 卡尔曼滤波不要求保存过去的测量数据，当新的数据到来时，根据新的数据和前一时刻的储值的估计，借助于系统本身的状态转移方程，按照一套递推公式，即可算出新的估值。卡尔曼递推算法大大减少了滤波装置的存储量和计算量，并且突破了平稳随机过程的限制，使卡尔曼滤波器适用于对时变信号的实时处理。

分形树__Matlab

%这是一个生成树的主函数，它的输入分别为每叉树枝的缩短比、树枝的偏角、生长次数. %注意：把这些程序全部保存在名为tree的M文件中再运行！！！！！！！！ %注意：把这些程序全部保存在名为tree的M文件中再运行！！！！！！！！ %注意：把这些程序全部保存在名为tree的M文件中再运行！！！！！！！！ %注意：把这些程序全部保存在名为tree的M文件中再运行！！！！！！！！ %注意：把这些程序全部保存在名为tree的M文件中再运行！！！！！！！！ %注意：把这些程序全部保存在名为tree的M文件中再运行！！！！！！！！ %注意：把这些程序全部保存在名为tree的M文件中再运行！！！！！！！！ %%小提示：若用做函数，请将虚线框内语句删去。 function f=tree(w,dtheata,NN) %%%--------------------虚线框--------------------%%% clear;clc;clf;w=0.8;dtheata=pi/6;NN=8;%建议生长次数NN不要超过10 %%%--------------------虚线框--------------------%%% n=2^NN;%从主枝算起，共需生成2^NN个树枝 for NNK=1:n x1=0; y1=0; r1=1; theata1=pi/2; dataway=ten2twoN(NNK,NN); %把每一个树枝的编号转化为一个NN位的二进制数 for NNL=1:NN if dataway(NNL)==0 [x2,y2,r2,theata2]=antmoveleft(x1,y1,r1,theata1,w,dtheata);%若路径数组上对应的数字为0，则向左生长 x1=x2; y1=y2; r1=r2; theata1=theata2; hold on %pause(eps) else [x2,y2,r2,theata2]=antmoveright(x1,y1,r1,theata1,w,dtheata);%否则，数字为1，向右生长 x1=x2; y1=y2; r1=r2; theata1=theata2; hold on %pause(eps) end end end hold off %--------------------------------------------------------------------------

基于matlab-的巴特沃斯低通滤波器的实现

基于matlab 的巴特沃斯低通滤波器的实现 一、课程设计的目的 运用MATLAB实现巴特沃斯低通滤波器的设计以及相应结果的显示，另外还对多种低通滤波窗口进行了比较。 二、课程设计的基本要求 1）熟悉和掌握MATLAB 的基本应用技巧。 2）学习和熟悉MATLAB相关函数的调用和应用。 3）学会运用MATLAB实现低通滤波器的设计并进行结果显示。 三、双线性变换实现巴特沃斯低通滤波器的技术指标： 1.采样频率10Hz。 2.通带截止频率fp=0.2*pi Hz。 3.阻带截止频率fs=0.3*pi Hz。 4.通带衰减小于1dB，阻带衰减大于20dB 四、使用双线性变换法由模拟滤波器原型设计数字滤波器 程序代码： T=0.1; FS=1/T; fp=0.2*pi;fs=0.3*pi; wp=fp/FS*2*pi; ws=fs/FS*2*pi; Rp = 1; % 通带衰减 As = 15; % 阻带衰减 OmegaP = (2/T)*tan(wp/2); % 频率预计 OmegaS = (2/T)*tan(ws/2); % 频率预计 %设计巴特沃斯低通滤波器原型

N = ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(OmegaP/OmegaS))); OmegaC = OmegaP/((10^(Rp/10)-1)^(1/(2*N))); [z,p,k] = buttap(N); %获取零极点参数 p = p * OmegaC ; k = k*OmegaC^N; B = real(poly(z)); b0 = k; cs = k*B; ds = real(poly(p)); [b,a] = bilinear(cs,ds,FS);% 双线性变换 figure(1);% 绘制结果 freqz(b,a,512,FS);%进行滤波验证 figure(2); % 绘制结果 f1=50; f2=250; n=0:63; x=sin(2*pi*f1*n)+sin(2*pi*f2*n); subplot(2,2,1);stem(x,'.'); title ('输入信号'); y=filter(b,a,x); subplot(2,2,2);stem(y,'.') ; title('滤波之后的信号'); figure(3) ; stem(y,'.') title（'输出的信号')）

卡尔曼滤波器及其简matlab仿真.

卡尔曼滤波器及其简matlab仿真 一、卡尔曼滤波的起源 谈到信号的分析与处理，就离不开滤波两个字。通常，信号的频谱处于有限的频率范围内，而噪声的频谱则散布在很广的频率范围内，为了消除噪声，可以进行频域滤波。但在许多应用场合，需要直接进行时域滤波，从带噪声的信号中提取有用信号。虽然这样的过程其实也算是对信号的滤波，但其所依据的理论，即针对随机信号的估计理论，是自成体系的。人们对于随机信号干扰下的有用信号不能“确知”，只能“估计”。为了“估计”，要事先确定某种准则以评定估计的好坏程度。 1960年卡尔曼发表了用递归方法解决离散数据线性滤波问题的论文A New Approach to Linear Filtering and Prediction Problems（线性滤波与预测问题的新方法），在这篇文章里一种克服了维纳滤波缺点的新方法被提出来，这就是我们今天称之为卡尔曼滤波的方法。卡尔曼滤波应用广泛且功能强大，它可以估计信号的过去和当前状态甚至能估计将来的状态即使并不知道模型的确切性质。 其基本思想是以最小均方误差为最佳估计准则，采用信号与噪声的状态空间模型利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计，求出当前时刻的估计值。算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的估计。 对于解决很大部分的问题，它是最优，效率最高甚至是最有用的。它的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 卡尔曼滤波不要求保存过去的测量数据，当新的数据到来时，根据新的数据和前一时刻的储值的估计，借助于系统本身的状态转移方程，按照一套递推公式，即可算出新的估值。卡尔曼递推算法大大减少了滤波装置的存储量和计算量，并且突破了平稳随机过程的限制，使卡尔曼滤波器适用于对时变信号的实时处理。 二、卡尔曼滤波的原理

FIR低通滤波器+matlab编程+滤波前后图形

Matlab实现振动信号低通滤波 附件txt中的数字是一个实测振动信号，采样频率为5000Hz，试设计一个长度为M=32的FIR低通滤波器，截止频率为600Hz，用此滤波器对此信号进行滤波。要求： （1）计算数字截止频率； （2）给出滤波器系数； （3）绘出原信号波形； （4）绘出滤波后的信号波形； 解答过程： 第一部分：数字截止频率的计算 =600/5000/2=0.24 数字截止频率等于截止频率除以采样频率的一半，即 n 第二部分：滤波器系数的确定 在matlab中输入如下程序，即可得到滤波器系数： n=32 Wn=0.24 b=fir1(n,Wn) 得到的滤波器系数b为 Columns 1 through 9 -0.0008 -0.0018 -0.0024 -0.0014 0.0021 0.0075 0.0110 0.0077 -0.0054 Columns 10 through 18 -0.0242 -0.0374 -0.0299 0.0087 0.0756 0.1537 0.2166 0.2407 0.2166 Columns 19 through 27 0.1537 0.0756 0.0087 -0.0299 -0.0374 -0.0242 -0.0054 0.0077 0.0110 Columns 28 through 33 0.0075 0.0021 -0.0014 -0.0024 -0.0018 -0.0008 第三部分：原信号波形 将附件4中的dat文件利用识别软件读取其中的数据，共1024个点，存在TXT 文档中，取名bv.txt，并复制到matlab的work文件夹。 在matlab中编写如下程序： x0=load('zhendong.txt'); %找到信号数据地址并加载数据。 t=0:1/5000:1023/5000; %将数据的1024个点对应时间加载

(完整word版)扩展卡尔曼滤波算法的matlab程序

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;

Newton分形的原理及Matlab实现

龙源期刊网 https://www.sodocs.net/doc/034380101.html, Ｎｅｗｔｏｎ分形的原理及Ｍａｔｌａｂ实现 作者：张健徐聪全付勇智 来源：《电脑知识与技术》2009年第24期 摘要:详细推导了复平面上Newton迭代法的原理和计算公式,用MATLAB编制程序实现了Newton迭代算法,得到了一些奇异、绚丽的分形图形。对《数学实验》课程有一定的参考价值。 关键词:Newton迭代法;分形;Matlab;数学实验 中图分类号:TP312文献标识码:A文章编号:1009-3044(2009)24-6997-03 The Principles of Newton Fractal and it's Realization Using MATLAB ZHANG Jian, XU Cong-quan, FU Yong-zhi (Department of Basic Courses, Southwest Forestry College, Kunming 650224, China) Abstract: The Principles and formulas of Newton fractal was explained,fractal graphics of Newton iteration was created using Matlab. Key words: newton iteration; fractal; Matlab; mathematical experimental 分形是非线性科学的一个重要分支,应用于自然科学和社会科学的众多领域。其中,分形图形以其奇异、绚丽多彩的特点,广泛应用于纺织印染、广告设计、装潢设计、计算机美术教学 等领域[1]。 很多分形图形都是用迭代的方式实现的,Newton迭代法就是其中的一种。由Newton迭代 法产生的分形图形称为Newton分形[2]。很多文献都对Newton分形进行了介绍,但都没有详细的计算公式和算法说明,读者很难编制相应程序。本文详细介绍了复平面上Newton迭代法的原理和计算公式,设计了相应的实现算法,并用Matlab编制程序实现了Newton分形的绘制,生成了一些奇异、瑰丽的分形图形。

基于MATLAB的巴特沃斯滤波器

数字信号处理课程设计 2015年 6 月25 日

目录 一．设计目的： (3) 二．设计要求： (3) 三．设计内容： (4) 3.1选择巴特涡斯低通数据滤波器及双线性变换法的原因 (4) 3.2巴特沃思低通滤波器的基本原理 (4) 3.3双线性变换法原理 (5) 3.4数字滤波器设计流程图 (7) 3.5数字滤波器的设计步骤 (7) 四．用matlab实现巴特沃斯低通数字滤波器的仿真并分析 (9) 4.1巴特沃斯低通数字滤波器技术指标的设置 (9) 4.2用matlab实现巴特沃斯低通数字滤波器的仿真 (9) 4.3波形图分析： (12) 五.总结与体会 (13) 六．附录参考文献 (14) 2

一．设计目的： 该课程设计是测控技术与仪器专业的必修课，开设课程设计的目的使学生掌握数字信号处理的基本概念和基本理论，能够利用辅助工具进行FIR和IIR数字滤波器的设计，进行一维信号的频谱分析，并进行仿真验证。加强实践教学环节，加强学生独立分析、解决问题的能力，培养学生动手能力和解决实际问题的能力，实现宽口径教育。 (1)理解低通滤波器的过滤方法。 (2)进一步熟悉低通滤波器的基本应用。 (3)用仿真工具matlab软件对设计的滤波器进行软件和硬件仿真。 (6)将对仿真结果进行比较，从而检验滤波器滤波性能的准确性。 二．设计要求： 地震发生时，除了会产生地震波，还会由地层岩石在断裂、碰撞过程中所发生的震动产生次声波。它的频率大约在每秒十赫兹到二十赫兹之间（可以用11Hz和15Hz的两个信号的和进行仿真，幅度可以分别设定为1、2）。大气对次声波的吸收系数很小，因此它可以传播的很远，而且穿透性很强。通过监测次声波信号可以监测地震的发生、强度等信息，因为自然界中广泛存在着各种次声波，这就对地震产生的次声波产生了干扰（可以用白噪声模拟，方差为5），需要采取一定的处理方法，才能检测到该信号，要求设计检测方案；并处理方法给出具体的软件（可以以51系列单片机、STM32F407、TMS320F28335或TMS320F6745为例）。 假设地震次声波信号为x,输入x=sin(2*π*11*t)+2*sin(2*π*15*t)和伴有白噪声的合成信号，经过滤波器后滤除15Hz以上的分量，即只保留x=sin(2*π*11*t)+2*sin(2*π*15*t)的分量信号，来验证设计的滤波器是否达到了设计要求。 3

几个分形matlab实现

几个分形的matlab实现 摘要：给出几个分形的实例，并用matlab编程实现方便更好的理解分形，欣赏其带来的数学美感 关键字：Koch曲线实验图像 一、问题描述： 从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的两边代替，形成山丘形图形如下 图1 在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的两条边代替，再次形成新的图形如此迭代，形成Koch分形曲线。 二、算法分析： 考虑由直线段（2个点）产生第一个图形（5个点）的过程。图1中，设 1 P和 5 P分别为 原始直线段的两个端点，现需要在直线段的中间依次插入三个点 2 P， 3 P， 4 P。显然 2 P位 于线段三分之一处， 4 P位于线段三分之二处， 3 P点的位置可看成是由 4 P点以 2 P点为轴心，逆时针旋转600而得。旋转由正交矩阵 ?? ? ? ? ? ? ? - = ) 3 cos( ) 3 sin( ) 3 sin( ) 3 cos( π π π π A 实现。 算法根据初始数据（ 1 P和 5 P点的坐标），产生图1中5个结点的坐标。结点的坐标数组形成一个2 5?矩阵，矩阵的第一行为 1 P的坐标，第二行为 2 P的坐标……，第五行为 5 P的坐标。矩阵的第一列元素分别为5个结点的x坐标，第二列元素分别为5个结点的y坐标。 进一步考虑Koch曲线形成过程中结点数目的变化规律。设第k次迭代产生的结点数为k n，第1 + k次迭代产生的结点数为 1+ k n，则 k n和 1+ k n中间的递推关系为3 4 1 - = +k k n n。

三、实验程序及注释： p=[0 0;10 0]; %P为初始两个点的坐标,第一列为x坐标,第二列为y坐标 n=2; %n为结点数 A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)]; %旋转矩阵 for k=1:4 d=diff(p)/3; %diff计算相邻两个点的坐标之差,得到相邻两点确定的向量 %则d就计算出每个向量长度的三分之一,与题中将线段三等分对应 m=4*n-3; %迭代公式 q=p(1:n-1,:); %以原点为起点,前n-1个点的坐标为终点形成向量 p(5:4:m,:)=p(2:n,:); %迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标 p(2:4:m,:)=q+d; %用向量方法计算迭代后处于4k+2位置上的点的坐标 p(3:4:m,:)=q+d+d*A'; %用向量方法计算迭代后处于4k+3位置上的点的坐标 p(4:4:m,:)=q+2*d; %用向量方法计算迭代后处于4k位置上的点的坐标 n=m; %迭代后新的结点数目 end plot(p(:,1),p(:,2)) %绘出每相邻两个点的连线 axis([0 10 0 10]) 四、实验数据记录： 由第三部分的程序，可得到如下的Koch分形曲线： 图2 五、注记： １．参照实验方法，可绘制如下生成元的Koch 分形曲线：

matlab对卡尔曼滤波的仿真实现

MATLAB 对卡尔曼滤波器的仿真实现 刘丹，朱毅，刘冰 武汉理工大学信息工程学院，武汉（430070） E-mail ：liudan_ina@https://www.sodocs.net/doc/034380101.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()(?

几个分形的matlab实现资料

几个分形的matlab 实现 摘要：给出几个分形的实例，并用matlab 编程实现方便更好的理解分形，欣赏其带来的 数学美感 关键字：Koch 曲线 实验 图像 一、问题描述： 从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的两边代替，形成山丘形图形如下 图1 在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的两条边代替，再次形成新的图形如此迭代，形成Koch 分形曲线。 二、算法分析： 考虑由直线段（2个点）产生第一个图形（5个点）的过程。图1中，设1P 和5P 分别为原始直线段的两个端点，现需要在直线段的中间依次插入三个点2P ，3P ，4P 。显然2P 位于线段三分之一处，4P 位于线段三分之二处，3P 点的位置可看成是由4P 点以2P 点为轴心，逆时针旋转600 而得。旋转由正交矩阵 ?????? ? ?-=)3cos()3sin()3sin()3cos(ππππA 实现。 算法根据初始数据（1P 和5P 点的坐标），产生图1中5个结点的坐标。结点的坐标数组形成一个25?矩阵，矩阵的第一行为1P 的坐标，第二行为2P 的坐标……，第五行为5P 的坐标。矩阵的第一列元素分别为5个结点的x 坐标，第二列元素分别为5个结点的y 坐标。 进一步考虑Koch 曲线形成过程中结点数目的变化规律。设第k 次迭代产生的结点数为k n ，第1+k 次迭代产生的结点数为1+k n ，则k n 和1+k n 中间的递推关系为341-=+k k n n 。 三、实验程序及注释：

p=[0 0;10 0]; %P为初始两个点的坐标,第一列为x坐标,第二列为y坐标 n=2; %n为结点数 A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)]; %旋转矩阵 for k=1:4 d=diff(p)/3; %diff计算相邻两个点的坐标之差,得到相邻两点确定的向量 %则d就计算出每个向量长度的三分之一,与题中将线段三等分对应 m=4*n-3; %迭代公式 q=p(1:n-1,:); %以原点为起点,前n-1个点的坐标为终点形成向量 p(5:4:m,:)=p(2:n,:); %迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标 p(2:4:m,:)=q+d; %用向量方法计算迭代后处于4k+2位置上的点的坐标 p(3:4:m,:)=q+d+d*A'; %用向量方法计算迭代后处于4k+3位置上的点的坐标 p(4:4:m,:)=q+2*d; %用向量方法计算迭代后处于4k位置上的点的坐标 n=m; %迭代后新的结点数目 end plot(p(:,1),p(:,2)) %绘出每相邻两个点的连线 axis([0 10 0 10]) 四、实验数据记录： 由第三部分的程序，可得到如下的Koch分形曲线： 图2 五、注记： １．参照实验方法，可绘制如下生成元的Koch 分形曲线： 图3