这个卡尔曼滤波程序哪位大哥可以帮我解释一下

这个卡尔曼滤波程序哪位大哥可以帮我解释一下,第1张

clear

N=200;%取200个数

w(1)=0;

w=randn(1,N);%产生一个1×N的行向量,第一个数为0,w为过程噪声(其和后边的v在卡尔曼理论里均为高斯白噪声)

x(1)=0;%状态x初始值

a=1;%a为状态转移阵,此程序简单起见取1

for k=2:N

x(k)=ax(k-1)+w(k-1); %系统状态方程,k时刻的状态等于k-1时刻状态乘以状态转移阵加噪声(此处忽略了系统的控制量)

end

V=randn(1,N);%测量噪声

q1=std(V);

Rvv=q1^2;

q2=std(x);

Rxx=q2^2; %此方程未用到Rxx

q3=std(w);

Rww=q3^2; %Rvv、Rww分别为过程噪声和测量噪声的协方差(此方程只取一组数方差与协方差相同)

c=02;

Y=cx+V;%量测方差,c为量测矩阵,同a简化取为一个数

p(1)=0;%初始最优化估计协方差

s(1)=0;%s(1)表示为初始最优化估计

for t=2:N

p1(t)=a^2p(t-1)+Rww;%p1为一步估计的协方差,此式从t-1时刻最优化估计s的协方差得到t-1时刻到t时刻一步估计的协方差

b(t)=cp1(t)/(c^2p1(t)+Rvv);%b为卡尔曼增益,其意义表示为状态误差的协方差与量测误差的协方差之比(个人见解)

s(t)=as(t-1)+b(t)(Y(t)-acs(t-1));%Y(t)-acs(t-1)称之为新息,是观测值与一步估计得到的观测值之差,此式由上一时刻状态的最优化估计s(t-1)得到当前时刻的最优化估计s(t)

p(t)=p1(t)-cb(t)p1(t);%此式由一步估计的协方差得到此时刻最优化估计的协方差

end

t=1:N;

plot(t,s,'r',t,Y,'g',t,x,'b');%作图,红色为卡尔曼滤波,绿色为量测,蓝色为状态

%整体来说,此卡尔曼程序就是一个循环迭代的过程,给出初始的状态x和协方差p,得到下一时刻的x和p,循环带入可得到一系列的最优的状态估计值,此方法通常用于目标跟踪和定位。

%本人研究方向与此有关,有兴趣可以交流下。

clc;clear all;

%归一化模拟切比雪夫I型低通滤波器的设计

Wp=2pi1000;Ws=2pi1500;rp=3;rs=30;%设计滤波器的参数    

wp=1;ws=Ws/Wp;                       %频带变换得到归一化滤波器

[N,wc]=cheb1ord(wp,ws,rp,rs,'s');    %计算滤波器阶数和3dB截止频率  

[z,p,k]=cheb1ap(N,rp);               %计算归一化滤波器的零极点

[b,a]=zp2tf(z,p,k);                  %计算归一化滤波器系统函数的系数

w0=0:005pi:2pi;

[h0,w0]=freqs(b,a,w0);               %求频率响应 

figure(1)

plot(w0,20log10(abs(h0)),'k');

title('归一化模拟且比雪夫I型低通滤波器');

xlabel('频率f/Hz');ylabel('幅度/dB');grid;

%一般低通切比雪夫低通滤波器的设计

[B,A]=lp2lp(b,a,Wp);                 %将归一化滤波器转换为一般模拟滤波器

w1=0:2pi100:2pi30000;

[h1,w1]=freqs(B,A,w1);

figure(2)

plot(w1/(2pi),20log10(abs(h1)),'k');

title('一般模拟且比雪夫I型低通滤波器');

xlabel('频率f/Hz');ylabel('幅度/dB');grid;

%冲激响应不变法设计低通巴特沃斯数字滤波器

Fs=10000;                            %采样频率 

Wp1=Wp/Fs;Ws1=Ws/Fs;

rp1=3;rs1=30;                        %设计滤波器的参数  

[N1,Wn]=cheb1ord(Wp,Ws,rp1,rs1,'s')  %计算滤波器阶数

[b1,a1]=cheby1(N1,rp1,Wn,'s');       %样本AF的系数函数的分子分母系数

[bz,az]=impinvar(b1,a1,Fs);          %冲激响应不变法从AF到DF变换

[C1,B1,A1]=dir2par(bz,az)            %直接形式转换为并联型

w2=[Wp1,Ws1];                        %数字临界频率

[H,f]=freqz(bz,az);                  %绘制数字滤波器的幅频特性和相频特性

[db,mag,pha,grd,f]=freqz_m(bz,az);   %扩展函数检验滤波器的其他特性

figure(3)

plot(f/pi,db,'k');

title('冲激响应不变法设计低通切比雪夫I型数字滤波器');

xlabel('角频率w/pi');ylabel('振幅/dB');

axis([0,035,-30,5]);grid;

%用设计好的滤波器对信号进滤波处理

figure(4)

f1=500;f2=4000;                      %输入信号的频率

N=100;                               %数据长度

dt=1/Fs;n=0:N-1;t=ndt;              %采样间隔和时间序列

x=sin(2pif1t)+05cos(2pif2t); %输入信号

subplot(2,1,1),plot(t,x),title('输入信号');

y=filtfilt(bz,az,x);                 %用滤波器进行滤波处理

y1=filter(bz,az,x);                  %进行滤波处理

subplot(2,1,2),plot(t,y,t,y1,':'),title('输出信号');xlabel('时间/s');

legend('filtfilt','filter')          %加图例

freqz_mm文件

function[db,mag,pha,grd,w]=freqz_m(b,a)

[H,w]=freqz(b,a,1000,'whole');

mag=abs(H);

db=20log10((mag+eps)/max(mag));

pha=angle(H);

grd=grpdelay(pha);

dir2parm文件

function [C,B,A] = dir2par(b,a)

M=length(b);N=length(a);

[r1,p1,C]=residuez(b,a);

p=cplxpair(p1,10000000eps);

I=cplxcomp(p1,p);

r=r1(I);

K=floor(N/2);B=zeros(K,2);A=zeros(K,3);

if K2==N;

    for i=1:2:N-2

        Brow=r(i:1:i+1,:);

        Arow=p(i:1:i+1,:);

        [Brow,Arow]=residuez(Brow,Arow,[]);

        B(fix((i+1)/2),:)=real(Brow);

        A(fix((i+1)/2),:)=real(Arow);

    end

    [Brow,Arow]=residuez(r(N-1),p(N-1),[]);

    B(K,:)=[real(Brow) 0];A(K,:)=[real(Arow) 0];

else

    for i=1:2:N-1

        Brow=r(i:1:i+1,:);

        Arow=p(i:1:i+1,:);

        [Brow,Arow]=residuez(Brow,Arow,[]);

        B(fix((i+1)/2),:)=real(Brow);

        A(fix((i+1)/2),:)=real(Arow);

    end

end

cplxcompm文件

function I=cplxcomp(p1,p2)

I=[];

for j=1:1:length(p2)

    for i=1:1:length(p1)

        if(abs(p1(i)-p2(j))<00001)

            I=[I,i];

        end

    end

end

I=I';

老师给出的程序,大部分,都不用动,

仅仅改动几个地址,即可符合题目要求。

把程序中的:

70H,改为 60H;

71H,改为 61H;

81H,改为 71H;

82H,改为 72H。

即可。

I=imread('lenabmp');

inf=imfinfo('lenabmp')

figure,imshow(I)

X=grayslice(I,64);

imshow(X,pink(64))

load trees

figure,image(10,10,X)

imwrite(X,map,'treesbmp');

imfinfo('treesbmp')

figure,imshow(X,map)

BW=im2bw(X,map,06);

figure,imshow(BW)

I=imread('lenabmp');

inf=imfinfo('lenabmp')

figure,imshow(I)

X=grayslice(I,64);

figure,imshow(X,pink(64))

A=imread('lenabmp');

imshow(A)

B=fftshift(fft2(A));

figure;

imshow(log(abs(B)),[8,10])

clc;

clear all;

I=imread('lenatif');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %用中值滤波,多维滤波,使用中心为-4,-8的拉普

% %拉斯滤波器,高斯低通滤波,拉普拉斯滤波器进行滤波处理

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure;%figure1

subplot(2,2,1);

imshow(I);

title('原始图像');

J=imnoise(I,'salt & pepper',004);%加椒盐噪声

title('加椒盐噪声');

subplot(2,2,2);

imshow(J);

K=medfilt2(J,[4,4])%进行中值滤波;

subplot(2,2,3);

imshow(K);

title('进行中值滤波');

h=ones(3,3)/9;%多维滤波

w=h;

g=imfilter(I,w,'conv','replicate')

subplot(2,2,4);

imshow(g);

title('多维滤波');

%使用中心为-4,-8的拉普拉斯滤波器,

w4=[1 1 1;

1 -4 1;

1 1 1];

w8=[1 1 1;

1 -8 1;

1 1 1];

f=im2double(I);

g4=f-imfilter(f,w4,'replicate');

g8=f-imfilter(f,w8,'replicate');

imshow(f);

figure;%figure2

subplot(1,2,1);

imshow(g4);

title('中心为-4的拉普拉斯滤波');

subplot(1,2,2);

imshow(g8);

title('中心为-8的拉普拉斯滤波');

h3=fspecial('gaussian',[3,3],05);%高斯低通滤波

figure;%figure3

B4=filter2(h3,I);

subplot(1,2,1);

imshow(B4,[ ]);

title('高斯低通滤波');

h4=fspecial('laplacian',0);%使用拉普拉斯滤波器

B5=filter2(h4,I);

subplot(1,2,2);

imshow(B5,[ ]);

title('拉普拉斯滤波器');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %从空域的角度进行亮度变换

% %把灰度等级是10-100的变化到10-255

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure;%figure4

subplot(2,2,1);

imshow(I);

title('原始图像');

J2=imadjust(I,[],[],05);% 增强低灰度级

subplot(2,2,2);

imshow(J2);

title('增强低灰度级');

J3=imadjust(I,[ ],[ ],2);%增强高灰度级

subplot(2,2,3);

imshow(J3);

title('增强高灰度级');

a1=100/255;%把灰度等级是10-100的变化到10-255

a2=255/255;

a3=10/255;

J2=imadjust(I,[a3,a1],[a3,a2],[]);

subplot(2,2,4);

imshow(J2);

title('把灰度等级是10-100的变化到10-255');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %从频域的角度进行亮度变换

% %fft2

% %由于能量主要集中在低频部分

% %所以对低频进行处理可以得到理想的效果

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I=imread('lenatif');

up=05;%设置处理频率上限

down=009%%设置处理频率下限

figure;%figure5

subplot(421);

imshow(I);

title('原始图像');

f=double(I);

subplot(4,2,3);

imshow(log(abs(f)),[]);

title('unit8转化为double');

g=fft2(f);

g=fftshift(g);

subplot(4,2,5);

imshow(log(abs(g)),[]);

title('FFT2变化后的图像');

[M,N]=size(g);% 转换数据矩阵

y1=max(max(abs(g)));%求出最大频率

y2=min(min(abs(g)));%%求出最小频率

y3=(y1-y2)up+y2;%设置滤波上限

y4=(y1-y2)down+y2;%%设置滤波下限

for i=1:M

for j=1:N

if (abs(g(i,j))<y4)

g(i,j)=g(i,j)^11;%对低频部分进行灰度增强

end

end

end

result=ifftshift(g);

J2=ifft2(result);

J3=uint8(abs(J2));

subplot(427);

imshow(J3,[ ]);

title('频域处理后的图像');

subplot(422)

imhist(I,64);

subplot(424)

imhist(f,64);

subplot(426)

imhist(g,64);

subplot(428)

imhist(J3,64);

以上就是关于这个卡尔曼滤波程序哪位大哥可以帮我解释一下全部的内容,包括:这个卡尔曼滤波程序哪位大哥可以帮我解释一下、帮忙编一个matlab的低通滤波器的程序、用汇编语言写去极值平均滤波程序等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!

欢迎分享,转载请注明来源:内存溢出

原文地址: https://outofmemory.cn/zz/9664752.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-04-30
下一篇 2023-04-30

发表评论

登录后才能评论

评论列表(0条)

保存