x = 0; % R = input('请输入过程噪声方差R的值: ');; % 测量噪声协方差
Q = input('请输入观测噪声方差Q的值: '); % 过程状态协方差
tf = 100; % 模拟长度
N = 100; % 粒子滤波器中的粒子个数
xhat = x;
P = 2;
xhatPart = x;
for i = 1 : N
u = rand;
qtempsum = 0;
for j = 1 : N
qtempsum = qtempsum + q(j);
if qtempsum >= u
xpart(i) = xpartminus(j);
break;
end
for i = 1 : N
xpart(i) = x + sqrt(P) * randn;
end
xArr = [x];
yArr = [x^2 / 20 + sqrt(R) * randn];
xhatArr = [x];
PArr = [P];
xhatPartArr = [xhatPart];
close all;
for k = 1 : tf
% 模拟系统
x = 0.5 * x + 25 * x / (1 + x^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn;
y = x^2 / 20 + sqrt(R) * randn;
% 粒子滤波器
for i = 1 : N
xpartminus(i) = 0.5 * xpart(i) + 25 * xpart(i) / (1 + xpart(i)^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn;
ypart = xpartminus(i)^2 / 20;
vhat = y - ypart;% 观测和预测的差
q(i) = (1 / sqrt(R) / sqrt(2*pi)) * exp(-vhat^2 / 2 / R);
end
% 平均每一个估计的可能性
qsum = sum(q);
for i = 1 : N
q(i) = q(i) / qsum;%归一化权值
end
end
end
xhatPart = mean(xpart);
% 在图片中表示出数据
xArr = [xArr x];
yArr = [yArr y];
xhatArr = [xhatArr xhat];
PArr = [PArr P];
xhatPartArr = [xhatPartArr xhatPart];
end
t = 0 : tf;
figure;
plot(t, xArr, 'b.', t, xhatPartArr, 'k-');
set(gca,'FontSize',10); set(gcf,'Color','yellow'); xlabel('time step'); ylabel('state');
legend('真实值', 'PF估计值');
xhatRMS = sqrt((norm(xArr - xhatArr))^2 / tf); xhatPartRMS = sqrt((norm(xArr - xhatPartArr))^2 / tf); disp(['PF估计误差均方值 = ', num2str(xhatPartRMS)]);
2017 年秋季学期研究生课程考核 （读书报告、研究报告） 考核科目:雷达系统导论 学生所在（系）:电子与信息工程学院 学生所在学科:电子与同学工程 学生姓名: 学号: 学生类别: 考核结果阅卷人 第 1 页（共页）
几种非线性滤波算法的介绍与性能分析 作者姓名：学号： 专业院系：电信学院电子工程系 电子邮件： 摘要—非线性滤波算法在雷达目标跟踪中有着重要的应用，对雷达的跟踪性能有着至关重要的影响。好的滤波算法有利于目标航迹的建立及保持，能够得到较精确的目标位置，为发现目标后的后续工作提供可靠的数据依据。本文重点介绍了雷达数据处理中的几种非线性滤波算法：扩展卡尔曼滤波（EKF）、不敏卡尔曼滤波（UKF）、粒子滤波（PF），并且给出了一个利用这三种算法进行数据处理的一个实例，通过这个实例对比分析了这三种算法的性能以及优劣。 关键字—非线性滤波算法；扩展卡尔曼滤波；不敏卡尔曼滤波；粒子滤波； I.概述（一级表题格式） 在雷达对目标进行跟踪前要先对目标进行检测。对于满足检测条件的目标就需要进行跟踪，在跟踪的过程中可以利用新获得的数据完成对目标的进一步检测比如去除虚假目标等，同时利用跟踪获得数据可以进一步完成对目标动态特性的检测和识别。因此对目标进行准确的跟踪是雷达性能的一个重要指标。在检测到满足条件的目标后，根据目标运动状态建立目标运动模型，然后对目标跟踪算法进行设计，这是雷达目标跟踪中的核心部分。 目前主要的跟踪算法包括线性自回归滤波，两点外推滤波、维纳滤波、- αβ滤波、加权最小二乘滤波、维纳滤波和卡尔曼滤波[1]。对于线性系统而言最优滤波的方法就是卡尔曼滤波，卡尔曼滤波是线性高斯模型下的最优状态估计算法。但是实际问题中目标的运动模型往往不是线性的，因此卡尔曼滤波具有很大的局限性。目前主要用的非线性滤波算法可以分为高斯滤波和粒子滤波[2]。不敏卡尔曼滤波和扩展卡尔曼滤波就是高斯滤波中的典型代表，也是应用相对较为广泛的。粒子滤波的应用范围比高斯滤波的适用范围要广，对于系统状态非线性，观测模型非高斯等问题都有很好的适用性。本文具体分析阐述了扩展卡尔曼滤波算法，不敏卡尔曼滤波算法，粒子滤波算法，并且通过一个实例利用仿真的方法分析了这三种算法在滤波性能上的优劣，最后对这三种算法做了一定的总结。 我本科毕业设计题目为《基于历史数据的路径生成算法研究》，由于我是跨专业保研到电信学院，该课题所研究内容不属于雷达系统研究范围，是一种城市路网最快路径生成算法。 II.几种非线性滤波算法 A.扩展卡尔曼滤波 扩展卡尔曼滤波是将非线性系统转换为近似的线性系统的一种方法，其核心思想是围绕滤波值将非线性函数展开成泰勒级数并略去二阶及以上的项，得到一个近似的线性化模型，然后应用卡尔曼滤波完成状态估计。 扩展卡尔曼滤波状态空间模型： k k k w x f+ = + ) ( x 1 状态方程 k k k v x h+ =) ( z观测方程 其中(.) f和(.) h为非线性函数 在扩展卡尔曼滤波中，状态的预测以及观测值的预测由非线性函数计算得出，线性卡尔曼滤波中的状态转移矩阵A阵和观测矩阵H阵由f和h函数的雅克比矩阵代替。 对 (.) f和(.) h Taylor展开，只保留一次项有： ) ? ( ) ?( ) ( k k k k k x x A x f x f- + ≈ ) ? ( ) ?( ) ( k k k k k x x H x h x h- + ≈ 其中： k k x x k k dx df A ?= =为f对 1- k x求导的雅克比矩阵 k k x x k k dx dh H ?= =为h对 1- k x求导的雅克比矩阵 ) ?( ? 1-k k x f x=，于是可以得出： k k k k k k k w x A x f x A x+ - + ≈ + ) ? ) ?( ( 1 k k k k k k k v x H x h x H z+ - + ≈ + ) ? ) ?( ( 1 通过以上变换，将非线性问题线性化。接下来EKF 滤波过程同线性卡尔曼滤波相同，公式如下: )) | (?( ) |1 ( X?k k X f k k= + ) ( ) ( ) | ( ) ( ) |1 (P k Q k k k P k k k+ Φ' Φ = + )1 ( )1 ( ) |1 ( )1 ( )1 (S+ + + ' + + = +k R k H k k P k H k )1 ( )1 ( ) |1 ( )1 ( K1+ + ' + = +-k S k H k k P k
数字滤波器matlab的源代码 function lvbo(Ua,Ub,choise) %参考指令：lvbo(2*pi,10*pi,1/0/-1) U1=min(Ua,Ub); U2=max(Ua,Ub); Us=16*U2; T=2*pi/Us; T_sum=4*max(2*pi/Ua,2*pi/Ub); sum=T_sum/T; t=T:T:T_sum; x=sin(U1*t)+0.8*sin(U2*t); X=DFT(x); figure(1); subplot(221) U=Us/sum:Us/sum:Us; stem(U,abs(X));grid on axis([Us/sum,Us/2,0,1.2*max(abs(X))]) title('原模拟信号采样频谱图') Ucd=U1+(U2-U1)*1/5;Usd=U2-(U2-U1)*1/5; switch choise case 1 Hz_ejw=IIR_DF_BW(Ucd,1,Usd,30,T,sum); case -1 Hz_ejw=IIR_DF_CF(Ucd,1,Usd,30,T,sum); case 0 Hz_ejw=FIR_DF_HM(U1,U2,T,sum); otherwise Hz_ejw=IIR_DF_BW(Ucd,1,Usd,30,T,sum); end Y=X.*Hz_ejw; y=1/sum*conj(DFT(conj(Y))); figure(1); subplot(224) plot(t,real(y)); title('模拟信号滤波后');grid on axis([0,T_sum,-max(real(y))*1.5,max(real(y))*1.5]) subplot(222); plot(t,x); hold on
毕 业 设 计 （论文） 开 题 报 告 姓名: 学号： 学院： 专业: 课题：基于粒子滤波的移动目标跟踪导师： 时间：
1.本课题研究的目的及意义： 粒子滤波(Particle Filter, PF)是一种基于蒙特卡罗(Monte Carlo)仿真的方法，它利用状态空间的一组带权值的随机样本(粒子)逼近状态变量的概率密度函数，每个样本代表系统的一个可能状态，可以得到状态的最小方差估计。粒子滤波算法摆脱了解决非线性滤波问题时随机量必须满足高斯分布的制约条件，因此，近几年来它在计算机视觉、目标跟踪、机器学习等领域受到了广泛的关注。另外，粒子滤波器的多模态处理能力，也是它应用广泛的原因之一。 本课题主要关注粒子滤波算法在目标跟踪领域的应用，随着计算机技术的发展，人们开始利用计算机来处理数字图像，包括图像增强，图像恢复，图像检索等等，而视频中运动目标的跟踪一直是计算机视觉、图像处理和模式识别等领域非常重要的研究课题。但是传统的目标跟踪方法存在着很多的局限性与不足之处，比如对非刚性目标跟踪时如何准确提取合适的目标特征进行跟踪，以及如何应对跟踪过程中的遮挡问题和复杂背景等等，也就难以保证跟踪的实时性和有效性。然而诸如此类的问题现在可以借助基于粒子仿真的方法来解决，在动态系统的模型选择，故障检测、诊断方面，出现了基于粒子的假设检验、粒子多模型、粒子似然度比检测等方法。同时，粒子滤波较之卡尔曼滤波（Kalman Filter）等在非线性非高斯系统领域中存在的优势，也决定了它的应用范围更加宽泛。 本课题旨在通过研究深入理解粒子滤波的原理及其算法，并利用MATLAB软件的图像处理功能，成功将粒子滤波算法应用于目标跟踪领域，最终实现对视频中运动目标的准确跟踪与检测。 2. 本课题国内外同类研究现状： 基于粒子滤波极强的实用性，国内外学者对此已经进行了大量研究，提出了许多用于跟踪的有效算法。这些方法主要可以分为两类：(1) 基于运动的方法：依据某种强健的算法，把一段时间内的具有运动一致性的点归为一类，如光流法和特征点法，但是计算量较大。(2) 基于模型的方法：主要依据高层的语义表示和知识描述来完成目标的跟踪。利用目标中信息部分的不同，可分为基于目标边界、基于目标区域的方法。但由于目标本身的信息较多，如不加简化，将不可避免地带来信息匹配时的大量运算。因此，对于实时性要求很高的运动目标的跟踪技术而言，如何选取目标的特征信息，并在可靠的前提下简化运算是目标跟踪的关键。本研究将在借鉴前人研究成果的基础
void CFastTrackingDlg::OnBnClickedButton4() { // TODO: Add your control notification handler code here //1.condensation setup const int stateNum=4; const int measureNum=2; const int sampleNum=2000; CvConDensation* condens = cvCreateConDensation(stateNum,measureNum,sampleNum); CvMat* lowerBound; CvMat* upperBound; lowerBound = cvCreateMat(stateNum, 1, CV_32F); upperBound = cvCreateMat(stateNum, 1, CV_32F); cvmSet(lowerBound,0,0,0.0 ); cvmSet(upperBound,0,0,winWidth ); cvmSet(lowerBound,1,0,0.0 ); cvmSet(upperBound,1,0,winHeight ); cvmSet(lowerBound,2,0,0.0 ); cvmSet(upperBound,2,0,0.0 ); cvmSet(lowerBound,3,0,0.0 ); cvmSet(upperBound,3,0,0.0 ); float A[stateNum][stateNum] ={ 1,0,1,0, 0,1,0,1, 0,0,1,0, 0,0,0,1 }; memcpy(condens->DynamMatr,A,sizeof(A)); cvConDensInitSampleSet(condens, lowerBound, upperBound); CvRNG rng_state = cvRNG(0xffffffff); for(int i=0; i < sampleNum; i++){ condens->flSamples[i][0] = float(cvRandInt( &rng_state ) % winWidth); //width condens->flSamples[i][1] = float(cvRandInt( &rng_state ) % winHeight);//height } CvFont font; cvInitFont(&font,CV_FONT_HERSHEY_SCRIPT_COMPLEX,1,1); char* winName="condensation"; cvNamedWindow(winName); cvSetMouseCallback(winName,mouseEvent); IplImage* img=cvCreateImage(cvSize(winWidth,winHeight),8,3); bool isPredictOnly=false;//trigger for prediction only,press SPACEBAR
m a t l a b程序之——滤波器(带通-带阻)
matlab程序之——滤波器（带通，带阻） 以下两个滤波器都是切比雪夫I型数字滤波器，不是巴特沃尔滤波器，请使用者注意！ 1.带通滤波器 function y=bandp(x,f1,f3,fsl,fsh,rp,rs,Fs) %带通滤波 %使用注意事项：通带或阻带的截止频率与采样率的选取范围是不能超过采样率的一半 %即，f1,f3,fs1,fsh,的值小于 Fs/2 %x:需要带通滤波的序列 % f 1：通带左边界 % f 3：通带右边界 % fs1：衰减截止左边界 % fsh：衰变截止右边界 %rp：边带区衰减DB数设置 %rs：截止区衰减DB数设置 %FS：序列x的采样频率 % f1=300;f3=500;%通带截止频率上下限 % fsl=200;fsh=600;%阻带截止频率上下限 % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值 % Fs=2000;%采样率 % wp1=2*pi*f1/Fs; wp3=2*pi*f3/Fs; wsl=2*pi*fsl/Fs; wsh=2*pi*fsh/Fs; wp=[wp1 wp3]; ws=[wsl wsh]; % % 设计切比雪夫滤波器； [n,wn]=cheb1ord(ws/pi,wp/pi,rp,rs); [bz1,az1]=cheby1(n,rp,wp/pi); %查看设计滤波器的曲线 [h,w]=freqz(bz1,az1,256,Fs); h=20*log10(abs(h));
figure;plot(w,h);title('所设计滤波器的通带曲线');grid on; y=filter(bz1,az1,x); end 带通滤波器使用例子 %-------------- %带通滤波器测试程序 fs=2000; t=(1:fs)/fs; ff1=100; ff2=400; ff3=700; x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t)+sin(2*pi*ff3*t); figure; subplot(211);plot(t,x); subplot(212);hua_fft(x,fs,1); % y=filter(bz1,az1,x); y=bandp(x,300,500,200,600,0.1,30,fs); figure; subplot(211);plot(t,y); subplot(212);hua_fft(y,fs,1); %调用到的hua_fft()函数代码如下 function hua_fft(y,fs,style,varargin) %当style=1,画幅值谱；当style=2,画功率谱;当style=其他的，那么花幅值谱和功率谱 %当style=1时，还可以多输入2个可选参数 %可选输入参数是用来控制需要查看的频率段的 %第一个是需要查看的频率段起点 %第二个是需要查看的频率段的终点 %其他style不具备可选输入参数，如果输入发生位置错误 nfft= 2^nextpow2(length(y));%找出大于y的个数的最大的2的指数值（自动进算最佳FFT步长nfft） %nfft=1024;%人为设置FFT的步长nfft y=y-mean(y);%去除直流分量 y_ft=fft(y,nfft);%对y信号进行DFT，得到频率的幅值分布 y_p=y_ft.*conj(y_ft)/nfft;%conj()函数是求y函数的共轭复数，实数的共轭复数是他本身。
完整的维纳滤波器Matlab源程序学术2009-11-20 20:31:58 阅读308 评论0 字号：大中小 完整的维纳滤波器Matlab源程序 clear;clc; %输入信号 A=1; %信号的幅值 f=1000; %信号的频率 fs=10^5; %采样频率 t=(0:999); %采样点 Mlag=100; %相关函数长度变量 x=A*cos(2*pi*f*t/fs); %输入正弦波信号 xmean=mean(x); %正弦波信号均值 xvar=var(x,1); %正弦波信号方差 xn=awgn(x,5); %给正弦波信号加入信噪比为20dB的高斯白噪声 figure(1) plot(t,xn) %绘制输入信号图像 title('输入信号图像') xlabel('x轴单位：t/s','color','b') ylabel('y轴单位：f/HZ','color','b') xnmean=mean(xn) %计算输入信号均值xnms=mean(xn.^2) %计算输入信号均方值xnvar=var(xn,1) %计算输入信号方差 Rxn=xcorr(xn,Mlag,'biased'); %计算输入信号自相关函数figure(2) subplot(221) plot((-Mlag:Mlag),Rxn) %绘制自相关函数图像 title('输入信号自相关函数图像') [f,xi]=ksdensity(xn); %计算输入信号的概率密度，f 为样本点xi处的概率密度 subplot(222) plot(xi,f) %绘制概率密度图像 title('输入信号概率密度图像') X=fft(xn); %计算输入信号序列的快速离散傅里叶变换 Px=X.*conj(X)/600; %计算信号频谱 subplot(223) semilogy(t,Px) %绘制在半对数坐标系下频
x = 0; % R = input('请输入过程噪声方差R的值: ');; % 测量噪声协方差 Q = input('请输入观测噪声方差Q的值: '); % 过程状态协方差 tf = 100; % 模拟长度 N = 100; % 粒子滤波器中的粒子个数 xhat = x; P = 2; xhatPart = x; for i = 1 : N u = rand; qtempsum = 0; for j = 1 : N qtempsum = qtempsum + q(j); if qtempsum >= u xpart(i) = xpartminus(j); break; end for i = 1 : N xpart(i) = x + sqrt(P) * randn; end xArr = [x]; yArr = [x^2 / 20 + sqrt(R) * randn]; xhatArr = [x]; PArr = [P]; xhatPartArr = [xhatPart]; close all; for k = 1 : tf % 模拟系统 x = 0.5 * x + 25 * x / (1 + x^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn; y = x^2 / 20 + sqrt(R) * randn; % 粒子滤波器 for i = 1 : N xpartminus(i) = 0.5 * xpart(i) + 25 * xpart(i) / (1 + xpart(i)^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn; ypart = xpartminus(i)^2 / 20; vhat = y - ypart;% 观测和预测的差 q(i) = (1 / sqrt(R) / sqrt(2*pi)) * exp(-vhat^2 / 2 / R); end % 平均每一个估计的可能性 qsum = sum(q); for i = 1 : N q(i) = q(i) / qsum;%归一化权值 end
数字信号处理
第一章概述 《数字信号处理》课程是通信专业的一门重要专业基础课，是信息的数字化处理、存储和应用的基础。通过该课程的课程设计实践，使我们对信号与信息的采集、处理、传输、显示、存储、分析和应用等有一个系统的掌握和理解，巩固和运用在《数字信号处理》课程中所学的理论知识和实验技能，掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。 数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。 其中，设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用
最广泛的是双线性变换法。 我们在课本中学到基本设计过程是： ①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标； ②设计过渡模拟滤波器； ③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。 而MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。第六章介绍的滤波器设计函数butter、cheby1 、cheby2 和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。 第二章总体方案设计 首先我将所给信号用MATLAB作图分析，然后通过观察st的幅频特性曲线，确定用高通滤波器作为处理信号的滤波器。选取滤波器的通带最大衰减为0.1dB,阻带最小衰减为60dB 为参数。 然后通过编程序调用MATLAB滤波器设计函数ellipord和ellip设计椭圆滤波器；通过编程序调用函数cheb1ord和cheby1设计切比雪夫滤波器，并绘图显示其幅频响应特性曲线。最后使用用滤波器实现函数filter，用两个滤波器分别对信号st进行滤波后绘图显示时域波形，观察滤波效果。 实验程序框图如图所示：
大庆石油学院学报 第32卷第3期2008年6月J OU RNAL OF DAQ IN G PETROL EUM INSTITU TE Vol.32No.3J un.2008 收稿日期:2007212224;审稿人:付光杰;编辑:郑丽芹 作者简介:任伟建(1963-),女,博士生导师,教授,主要从事复杂系统的控制及故障诊断方面的研究. 基于粒子滤波和贝叶斯估计的目标跟踪 任伟建1,山茂泉1,谢 锋2,王文东3 (1.大庆石油学院电气信息工程学院,黑龙江大庆 163318; 2.大庆油田有限责任公司第二采油厂,黑龙江大庆 163414; 3.大庆钻井技术服务公司钻井工具分公司,黑龙江大庆 163461) 摘 要:针对颜色直方图的彩色物体的运动目标,在各种噪声的干扰下多呈现非线性和非高斯的特点,利用粒子滤 波的方法进行运动估计和跟踪.利用粒子滤波对非线性和非高斯的有效逼近的性质,获得粒子的后验概率分布,估计目 标状态,实现目标的有效跟踪.采用累加权值概率并且引入随机正态分布进行采样,保证粒子的多样性,有效避免粒子退 化问题.仿真结果表明该方法的有效性. 关 键 词:粒子滤波;贝叶斯估计;目标跟踪;彩色直方图 中图分类号:TP182 文献标识码:A 文章编号:100021891(2008)0320067204 0 引言 目标存在变化多样和跟踪设备对环境适应性不完善等问题,复杂环境下的运动目标跟踪是个难题[1,2].为了有效跟踪运动目标,必须对运动对象进行有效的估计,利用已有的信息,获得当前运动物体估计状态,然后利用现有观察数据对运动状态进行修正.该类问题经常采用广义卡尔曼滤波方法.广义卡尔曼滤波依赖于模型的线性化和高斯假设.在估计系统状态和方差时,由于线性逼近,可能导致滤波发散.且如果密度函数不是高斯分布,该方法估计精度不高.近年来出现一种新的最优非线性方法———粒子滤波,它源自序列蒙特卡罗方法[3].该方法不受动态系统各个随机变量的限制,能够有效地应用于非线性、非高斯的运动系统中. 文中首先对选定区域目标建立颜色直方图模型,然后在选定区域附近产生目标粒子区域,利用巴特查理亚系数测量粒子区域和选定区域2种分布之间的相似度,运用粒子滤波估计方法实现运动目标的跟踪.在跟踪过程中,粒子存在退化现象.文献[4]采取重采样方法在一定程度上解决了退化问题,但由于重采样是根据权值大小进行的,导致采样后的粒子由大量重复的粒子构成,失去了多样性.文中采取概率累加的方法保持粒子的多样性,防止粒子退化,取得较好的效果. 1 运动目标模型 在确定运动目标后,建立基于指数分布的统计模型.在区域中心,属于运动目标的概率为1,在偏离中心的距离大于阈值时,概率属于指数衰减[5]: p pos (z i )=1,‖z i ‖≤T ; exp -‖z i ‖-T max (‖z i ‖-T )N i =1 ,‖z i ‖>T ,(1)可得到目标的统计直方图分布模型: p pos (u )=C 6N i =1p pos (z i )δ(b (z i )-u ),(2) C =1 6N i =1p pos (z i ). (3)
DSP 设计滤波器报告 姓名：张胜男 班级：07级电信（1）班 学号：078319120 一·低通滤波器的设计 （一）实验目的：掌握IIR 数字低通滤波器的设计方法。 （二）实验原理： 1、滤波器的分类 滤波器分两大类：经典滤波器和现代滤波器。 经典滤波器是假定输入信号)(n x 中的有用成分和希望取出的成分各自占有不同的频带。这样，当)(n x 通过一个线性系统（即滤波器）后可讲欲去除的成分有效的去除。 现代滤波器理论研究的主要内容是从含有噪声的数据记录（又称时间序列）中估计出信号的某些特征或信号本身。 经典滤波器分为低通、高通、带通、带阻滤波器。每一种又有模拟滤波器（AF ）和数字滤波器（DF ）。对数字滤波器，又有IIR 滤波器和FIR 滤波器。 IIR DF 的转移函数是： ∑∑=-=-+==N k k k M r r r z a z b z X z Y z H 10 1)()()( FIR DF 的转移函数是： ∑-=-=10)()(N n n z n h z H FIR 滤波器可以对给定的频率特性直接进行设计，而IIR 滤波器目前最通用的方法是利用已经很成熟的模拟滤波器的设计方法进行设计。 2、滤波器的技术要求 低通滤波器： p ω：通带截止频率（又称通带上限频率） s ω：阻带下限截止频率 p α：通带允许的最大衰减 s α：阻带允许的最小衰减 （p α，s α的单位dB ） p Ω：通带上限角频率 s Ω：阻带下限角频率 （s p p T ω=Ω，s s s T ω=Ω）即 C p p F ωπ2=Ω C s s F ωπ2=Ω 3、IIR 数字滤波器的设计步骤：
clear; % tic; x = 0.1; % 初始状态 x_estimate = 1;%状态的估计 e_x_estimate = x_estimate; %EKF的初始估计 u_x_estimate = x_estimate; %UKF的初始估计 p_x_estimate = x_estimate; %PF的初始估计 Q = 10;%input('请输入过程噪声方差Q的值: '); % 过程状态协方差 R = 1;%input('请输入测量噪声方差R的值: '); % 测量噪声协方差 P =5;%初始估计方差 e_P = P; %UKF方差 u_P = P;%UKF方差 pf_P = P;%PF方差 tf = 50; % 模拟长度 x_array = [x];%真实值数组 e_x_estimate_array = [e_x_estimate];%EKF最优估计值数组 u_x_estimate_array = [u_x_estimate];%UKF最优估计值数组 p_x_estimate_array = [p_x_estimate];%PF最优估计值数组 u_k = 1; %微调参数 u_symmetry_number = 4; % 对称的点的个数 u_total_number = 2 * u_symmetry_number + 1; %总的采样点的个数 linear = 0.5; N = 500; %粒子滤波的粒子数 close all; %粒子滤波初始N 个粒子 for i = 1 : N p_xpart(i) = p_x_estimate + sqrt(pf_P) * randn; end for k = 1 : tf % 模拟系统 x = linear * x + (25 * x / (1 + x^2)) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn; %状态值y = (x^2 / 20) + sqrt(R) * randn; %观测值 clear; % tic; x = 0.1; % 初始状态 x_estimate = 1;%状态的估计 e_x_estimate = x_estimate; %EKF的初始估计 u_x_estimate = x_estimate; %UKF的初始估计 p_x_estimate = x_estimate; %PF的初始估计 Q = 10;%input('请输入过程噪声方差Q的值: '); % 过程状态协方差
MATLAB大作业 院（系）：信息工程学院 专业：09通信工程 班级：通信一班 学生：钟锦慧 学号：20090610080118 指导教师：邹丹 2011年12月18日
MATLAB在滤波器设计中的应用 1. 绪论 从20世纪初至今，在通信与电子系统中，滤波器的研究和应用经历了漫长、艰辛而曲折的道路，滤波器在信号传输与信号处理中的重要地位和作用已经非常明显，所以滤波器的分析与设计更是应该重点研究的问题。滤波器，顾名思义，是对波进行过滤的器件。“波”是一个非常广泛的物理概念，在电子技术领域，“波”被狭义地局限于特指描述各种物理量的取值随时间起伏变化的过程。该过程通过各类传感器的作用，被转换为电压或电流的时间函数，称之为各种物理量的时间波形，或者称之为信号。因为自变量时间是连续取值的，所以称之为连续时间信号，又习惯地称之为模拟信号(Analog Signal)。随着数字式电子计算机(一般简称计算机)技术的产生和飞速发展，为了便于计算机对信号进行处理，产生了在抽样定理指导下将连续时间信号变换成离散时间信号的完整的理论和方法。也就是说，可以只用原模拟信号在一系列离散时间坐标点上的样本值表达原始信号而不丢失任何信息，波、波形、信号这些概念既然表达的是客观世界中各种物理量的变化，自然就是现代社会赖以生存的各种信息的载体。信息需要传播，靠的就是波形信号的传递。信号在它的产生、转换、传输的每一个环节都可能由于环境和干扰的存在而畸变，有时，甚至是在相当多的情况下，这种畸变还很严重，以致于信号及其所携带的信息被深深地埋在噪声当中了[。 2. MATLAB简介 2.1 MATLAB的概述 20世纪70年代，美国新墨西哥大学计算机科学系主任Cleve Moler为了减轻学生编程的负担，用FORTRAN编写了最早的MA TLAB。1984年由Little、Moler、Steve Bangert合作成立了的MathWorks公司正式把MA TLAB推向市场。到20世纪90年代，MA TLAB已成为国际控制界的标准计算软件。 MA TLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MA TLAB和Simulink两大部分。 MA TLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测
function ParticleEx1 % Particle filter example, adapted from Gordon, Salmond, and Smith paper. x = 0.1; % initial state Q = 1; % process noise covariance R = 1; % measurement noise covariance tf = 50; % simulation length N = 100; % number of particles in the particle filter xhat = x; P = 2; xhatPart = x; % Initialize the particle filter. for i = 1 : N xpart(i) = x + sqrt(P) * randn; end jArr = [0]; xArr = [x]; yArr = [x^2 / 20 + sqrt(R) * randn]; xhatArr = [x]; PArr = [P]; xhatPartArr = [xhatPart]; close all; for k = 1 : tf % System simulation x = 0.5 * x + 25 * x / (1 + x^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn;%状态方程 y = x^2 / 20 + sqrt(R) * randn;%观测方程 % Extended Kalman filter F = 0.5 + 25 * (1 - xhat^2) / (1 + xhat^2)^2; P = F * P * F' + Q; H = xhat / 10; K = P * H' * (H * P * H' + R)^(-1); xhat = 0.5 * xhat + 25 * xhat / (1 + xhat^2) + 8 * cos(1.2*(k-1));%预测 xhat = xhat + K * (y - xhat^2 / 20);%更新 P = (1 - K * H) * P; % Particle filter for i = 1 : N xpartminus(i) = 0.5 * xpart(i) + 25 * xpart(i) / (1 + xpart(i)^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn;
%高斯低通滤波器 RGB=imread('132.jpg'); I0=rgb2gray(RGB); subplot(2,3,1),imshow(I0);title('原图'); I1=imnoise(I0,'gaussian');%对原图像加噪声 subplot(2,3,2),imshow(I1);title('加入噪声后') %将灰度图像的二维不连续Fourier变换的零频率成分移到频谱的中心 s=fftshift(fft2(I1)); subplot(2,3,3),imshow(log(1+abs(s)),[]);title('fftshift'); [M,N]=size(s);%分别返回s的行数到M中，列数到N中 %GLPF滤波 d0=50;%初始化d0 n1=floor(M/2);%对M/2进行取整 n2=floor(N/2);%对N/2进行取整 for i=1:M for j=1:N d=sqrt((i-n1)^2+(j-n2)^2);%点（i,j）到傅立叶变换中心的距离 h(i,j)=1*exp(-1/2*(d^2/d0^2));%GLPF滤波函数
s(i,j)=h(i,j)*s(i,j);%GLPF滤波后的频域表示 end end s=ifftshift(s);%对s进行反FFT移动 %对s进行二维反离散的Fourier变换后，取复数的实部转化为无符号8位整数 s=uint8(real(ifft2(s))); subplot(2,3,4),imshow(h);title('传递函数');%显示GHPF滤波器的传递函数 subplot(2,3,5),imshow(s);title('GLPF滤波(d0=50)');%显示GLPF滤波处理后的图像
各类滤波器的MATLAB程序 一、理想低通滤波器 IA=imread('lena.bmp'); [f1,f2]=freqspace(size(IA),'meshgrid'); Hd=ones(size(IA)); r=sqrt(f1.^2+f2.^2); Hd(r>0.2)=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('lena.bmp'); [f1,f2]=freqspace(size(IA),'meshgrid'); Hd=ones(size(IA)); r=sqrt(f1.^2+f2.^2); Hd(r<0.2)=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('lena.bmp'); [f1,f2]=freqspace(size(IA),'meshgrid'); D=0.3; 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)); figure subplot(2,2,1),imshow(uint8(IA)); subplot(2,2,2),imshow(uint8(Ia)); figure surf(Hd,'Facecolor','interp','Edgecolor','none','Facelighting','phong');
虽然matlab提供了很多产生低通滤波器的函数，而且也提供了将低通转换为高通、带通等滤波器的方法函数，以及数字化的函数。但是为了简化设计及设计者方便考虑，matlab还提供了更为简便的产生各种滤波器的方法。 1 besself 功能：贝塞尔（Bessel）模拟滤波器设计。 格式：[b,a] = besself(n,Wn) [b,a] = besself(n,Wn,'ftype') [z,p,k] = besself(...) [A,B,C,D] = besself(...) 说明：besself函数可以设计模拟低通、高通、带通和带阻贝塞尔（Bessel）滤波器。 [b,a] = besself(n,Wn)返回截止频率为Wn（单位为弧度/秒）的n阶贝塞尔模拟 低通滤波器，b、a分别为滤波器传递函数的分子和分母系数向量（降幂排列）。 当Wn为二元向量，即Wn=[W1 W2] (W1 2 butter 功能：巴特沃思（Butterworth）模拟／数字滤波器设计。 格式：[b,a] = butter(n,Wn) [b,a] = butter(n,Wn,'ftype') [b,a] = butter(n,Wn,'s') [b,a] = butter(n,Wn,'ftype','s') [z,p,k] = butter(...) [A,B,C,D] = butter(...) 说明：butter函数可以设计模拟或数字的低通、高通、带通和带阻Butterworth 滤波 器。Butterworth滤波器可以使通带内的幅频响应最大程度地平坦，但这也使 得它的通带到阻带的过渡过程较慢。在这方面Chebyshev滤波器和椭圆滤波 器性能较好。 在设计数字滤波器时，butter函数中的参数Wn与besself函数有很大的区别， 它是一个相对量，其定义区间为Wn∈[0,1]，其中1对应于0.5fs，fs为采样 频率（单位Hz）；在设计模拟滤波器时，Wn采用真实频率，单位为Hz。 [b,a] = butter(n,Wn)返回截止频率为Wn（单位为弧度/秒）的n阶Butterworth 数字低通滤波器，b、a分别为滤波器传递函数的分子和分母系数向量（降幂 排列）。 当Wn为二元向量，即Wn = [W1 W2] (W1 % EKF UKF PF 的三个算法 clear; % tic; x = 0.1; % 初始状态 x_estimate = 1;%状态的估计 e_x_estimate = x_estimate; %EKF的初始估计 u_x_estimate = x_estimate; %UKF的初始估计 p_x_estimate = x_estimate; %PF的初始估计 Q = 10;%input('请输入过程噪声方差Q的值: '); % 过程状态协方差 R = 1;%input('请输入测量噪声方差R的值: '); % 测量噪声协方差 P =5;%初始估计方差 e_P = P; %UKF方差 u_P = P;%UKF方差 pf_P = P;%PF方差 tf = 50; % 模拟长度 x_array = [x];%真实值数组 e_x_estimate_array = [e_x_estimate];%EKF最优估计值数组 u_x_estimate_array = [u_x_estimate];%UKF最优估计值数组 p_x_estimate_array = [p_x_estimate];%PF最优估计值数组 u_k = 1; %微调参数 u_symmetry_number = 4; % 对称的点的个数 u_total_number = 2 * u_symmetry_number + 1; %总的采样点的个数 linear = 0.5; N = 500; %粒子滤波的粒子数 close all; %粒子滤波初始N 个粒子 for i = 1 : N p_xpart(i) = p_x_estimate + sqrt(pf_P) * randn; end for k = 1 : tf % 模拟系统 x = linear * x + (25 * x / (1 + x^2)) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn; %状态值 y = (x^2 / 20) + sqrt(R) * randn; %观测值 %扩展卡尔曼滤波器 %进行估计第一阶段的估计 e_x_estimate_1 = linear * e_x_estimate + 25 * e_x_estimate /(1+e_x_estimate^2) + 8 * cos(1.2*(k-1)); e_y_estimate = (e_x_estimate_1)^2/20; %这是根据k=1时估计值为1得到的观测值；只是这个由我估计得到的第24行的y也是观测值不过是由加了噪声的真实值得到的 %相关矩阵 e_A = linear + 25 * (1-e_x_estimate^2)/((1+e_x_estimate^2)^2);%传递矩阵 e_H = e_x_estimate_1/10; %观测矩阵 %估计的误差 e_p_estimate = e_A * e_P * e_A' + Q; %扩展卡尔曼增益 e_K = e_p_estimate * e_H'/(e_H * e_p_estimate * e_H' + R); %进行估计值的更新第二阶段 e_x_estimate_2 = e_x_estimate_1 + e_K * (y - e_y_estimate); %更新后的估计值的误差 e_p_estimate_update = e_p_estimate - e_K * e_H * e_p_estimate;EKF,UKF,PF三种算法的比较粒子滤波matlab仿真程序EKF,UKF,PF三种算法的比较matlab仿真程序
相关主题
文本预览