卡尔曼滤波用于语音增强的一个经典代码matlab-副本.docx

卡尔曼滤波用于语音增强的一个经典代码matlab -副本function output=KalmanSignalDenoiser(Noisy,Clean,fs)

%OUTPUT=KALMANSIGNALDENOISER(NOISY,CLEAN,FS)

%this purpose of this function is to demonstrate the capability of

kalman % filter for denoising noisy speech (corrupted by white noise).

Kalman % filtering of noisy speech usually have two steps:

%1 . Estimating the AR parameters of speech segment

%2 . Filtering the segment

%There are different approaches for extracting AR parameters of

noisy % speech in the literature, however, in this function none is

implemented. % Clean speech signal should be provided for this purpose.

%

%ARGUMENTS

%NOISY : noise contaminated speech

%CLEAN : clean speech signal

%FS : Sampling frequency which should be the same for both signals

%

%Output is the denoised speech signal

%

%Required functions:

%SEGMENT

%Sep-04

%Esfandiar Zavarehei

W=fix(.025*fs); %Window length is 25 ms

SP=1; %Shift percentage is 40% (10ms) %Overlap-Add method works good with this value(.4)

SpecP=13;

Window=ones(W,1);

x=segment(Clean,W,SP,Window);

y=segment(Noisy,W,SP,Window);

n=segment(Noisy-Clean,W,SP,Window);

R=var(n);

H=[zeros(1,SpecP-1) 1];

G=H'; GGT=G*H;

FUpper=[zeros(SpecP-1,1) eye(SpecP-1)];

I=eye(SpecP);

[A Q]=lpc(x,SpecP);

P=diag(repmat(R(1),1,SpecP));

o=zeros(1,W*size(x,2));% allocating memory to the output in advance

save a lot of

computation time

o(1:SpecP)=y(1:SpecP,1)';

hwb = waitbar(0,'Please wait...','Name','Processing'); start=SpecP+1;

Sp=o(1:SpecP)';

t=SpecP+1;

for n=1:size(x,2)

waitbar(n/size(x,2),hwb,['Please wait... '

num2str(fix(100*n/size(x,2)))

' %'])

F=[FUpper; fliplr(-A(n,2:end))];

for i=start:W

S_=F*Sp;

e=y(i,n)-S_(end);%innovation

P_=F*P*F'+GGT*Q(n);

K=(P_*H')/(H*P_*H' + R(n));

SOut=S_+K*e;

o(t-SpecP+1:t)=SOut'; %Notice that the previous SpecP-1 output samples are

updated again

P=(I-K*H)*P_;

Sp=SOut;

t=t+1;

end

start=1;

end

close(hwb)

output=o;

function FormantTrack=FTrackSig(signal,fs)

%F=FTRACKSIG(S)

%Formant Tracking

% Inputs: S, signal, the speech signal with Motorola Format

and sampling

%frequency fs (default=16000 KHz)

%Outputs: F, FormantTrack is a matrix each row of which is a

track of one of

%the formants

%

%Important Variables in the program: There are several

important variables

%in the program:

%Main Function (FTrackSig)

%MaxFormantFreq : The maximum frequency of the formant % MaxBWFreq : The maximum BandWidth of the formant % SegLength: Number of samplein each frame,

%SP: Shift Percentage of frames, (1-SP is overlap percentage) % NumberOfFormants: The number of formants to be chosen % PreEmphFact:

the pre-emphasis factor

%fs: sampling frequency

%FormantCand function

%P: order of LPC model

%----------------------------------

%What does this program do:

%The speech signal is first choped into segments and the pre-

emphasize, % then windowed. Then the LPC Model for these segments are

extracted. using % these coefficients the roor=ts of the LP model is calculated, sorted by % increasing frequency and the extreme poles

are eliminated. (an extreme % ploe is one with zero frequency or out of range frequency, or rather % large bandwidth). if the number of poles

are not sufficient the segments % will be re modeled using a larger

order for LPC (THIS IS STILL TO BE % DISCUSSED ACTUALLY). these sorted poles are then grouped with regards to % the number of poles required. if there were more than one group each % group will be evaluated

with regards to the different distributions of % phonems which are

loaded from another file, and then evaluated in terms % of the

trajectory continueity. the best formant will be chosen then using %

these evaluations.

%

%

%Author Esfandiar Zavarehei

%Date 5-Nov-2003

%Revised: 19-Nov-2003

%Revised: 21-Nov-2003

%

%

%Notes:

%The MAXIMUM FREQUENCY is a crucial parameter if you have a % low sampling frequency you have to adjust it your self. % The NUMBER OF

FORMANTS is 5 if you have a low sampling frequency % try in

decrease this number.

if nargin<2

fs=16000; % The Sampling frequency

end

SegLength=.025*fs; % This is equal to 25ms

SP=.4; % Shift Percentage

PreEmphFact=-0.98;

NumberOfFormants=5;

MaxFormantFreq=5000;

MaxBWFreq=600;

FreqRes=512;

sigLength=length(signal);

NumSeg=fix(sigLength-SegLength+SegLength*SP)/(SegLength*SP); signal=real(filter([1 PreEmphFact],1,signal)); %Pre-emphasizing BegPtr=1;

FCFreq=[];

FCBW=[];

hamwin=hamming(SegLength);

AddP=1;

for n=1:NumSeg

segment=signal(BegPtr:(BegPtr+SegLength-1)); %Chopping

segment=segment.*hamwin; %Windowing

segment=[segment;zeros((FreqRes-SegLength),1)]; %adjusting frequency

resoloution

Candidates=FormantCand(segment);

FCseg=angle(Candidates)*fs/(2*pi);

FCBW=abs(log(abs(Candidates))*fs/pi);

FCseg=FCseg(find(FCBW

FCseg=FCseg(find(FCseg

while (NumberOfFormants>length(FCseg))

Candidates=FormantCand(segment,AddP);

FCseg=angle(Candidates)*fs/(2*pi);

FCBW=abs(log(abs(Candidates))*fs/pi);

FCseg=FCseg(find(FCBW

FCseg=FCseg(find(FCseg

AddP=AddP+1;

end

if (n>1)

Formants=distFormant(FCseg,NumberOfFormants,LastFormant); else

Formants=distFormant(FCseg,NumberOfFormants);

end

ForMat=Formants.FM;

chance=1./Formants.dst;

[a b]=max(chance);

LastFormant=ForMat(:,b);

FormantTrack(:,n)=LastFormant;

AddP=1;

BegPtr=BegPtr+fix(SegLength*SP); %Moving the pointer to next

segment end

%-----------------------------extracting, and choosing the

appropriate poles of each segment

function Candidates=FormantCand(segment,AddP)

%5 NOV 2003

%This Program returns the Formant candidate of a segment %the segment, which is probably a piece of speech signal, %is not windowed, pre-emphasized, zero padded, etc. in this program. if (nargin<2) AddP=0; %The number to be added to the LP Order, if it did not have enough Poles

extracted

end

P=15+AddP; %number of LPC Coefitions

lpccof=lpc(segment,P);

rts=roots(lpccof);

Phase=angle(rts);

Amp=abs(rts);

posPhaseIndex=find((Phase>0.000001) & (Phase<(pi-0.000001))); Phase=Phase(posPhaseIndex);

Amp=Amp(posPhaseIndex);

[Phase sIndex]=sort(Phase);

Amp=Amp(sIndex);

%Candidates

Candidates=Amp.*exp(j*Phase);

%-------------------------------------------DISTFORMANT function

dist=distFormant(Poles,numberOfFormants,LastChosenFormant) %THIS

FUNCTION CALCULATE THE DISTANCE THE LAST FORMANT WITH ALL AVAILABLE %CANDIDATES AND RETURN A VECTOR (ROW VWCTOR) THAT CONTAINS THESE DISTANCES. %Formants are row vectors in this program

%if (length(Poles)

%dist.FM=[Poles;NaN*ones(numberOfFormants-length(Poles),1)];

%dist.dst=1;

%return

%end

FormantMat=(nchoosek(Poles,numberOfFormants))';

numberOfCands=size(FormantMat,2);

if (nargin==3)

dst=sum((FormantMat-LastChosenFormant(:,ones(1,numberOfCands))).^2);

elseif (nargin==2)

dst=ones(1,size(FormantMat,2));

else

error('not appropriate number of inputs') end

dist.FM=FormantMat;

dist.dst=dst;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function vec=mat2vec(mat)

%VEC=MAT2VEC(MAT)

%This function augments the rows of the input matrix MAT

%forming a vector VEC

%

%Example:

%

%A=[1216

%2419];

%B = mat2vec(MAT);

%Then Be will be:

%B=[12162419];

%

%Sep-04

%Esfandiar Zavarehei

numofrow=size(mat,1);

vec=[];

for i=1:numofrow

vec=[vec mat(i,:)];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function ReconstructedSignal=OverlapAdd(XNEW,yphase,windowLen,ShiftLen);

%Y=OverlapAdd(X,A,W,S);

%Y is the signal reconstructed signal from its spectrogram. X is

a matrix

%with each column being the fft of a segment of signal. A is the phase %angle of the spectrum which should have the same dimension as X. if it is %not given the phase angle of X is used which in the case of

real values is %zero (assuming that its the magnitude). W is the window length of time %domain segments if not given the length is assumed to be twice as long as %fft window length. S is the shift length of the segmentation process ( for %example in the case of non overlapping signals it is equal to W and in the %case of %50 overlap is equal to W/2. if not givven W/2 is used. Y is the %reconstructed time domain signal.

%Sep-04

%Esfandiar Zavarehei

if nargin<2

yphase=angle(XNEW);

end

if nargin<3

windowLen=size(XNEW,1)*2;

end

if nargin<4

ShiftLen=windowLen/2;

end

[FreqRes FrameNum]=size(XNEW);

sig=zeros((FrameNum-1)*ShiftLen+windowLen,1);

for i=1:FrameNum

start=(i-1)*ShiftLen+1;

spec=[XNEW(:,i).*exp(j*yphase(:,i))];

spec=[spec;flipud(conj(spec))];

sig(start:start+windowLen-1)=sig(start:start+windowLen-

1)+real(ifft(spec,windowLen));

end

ReconstructedSignal=sig;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555 function Seg=segment(signal,W,SP,Window)

%SEGMENT chops a signal to overlapping windowed segments

%A= SEGMENT(X,W,SP,WIN) returns a matrix which its columns are segmented % and windowed frames of the input one dimentional signal, X.

W is the % number of samples per window, default value W=256. SP is the

shift % percentage, default value SP=0.4. WIN is the window that is

multiplied by % each segment and its length should be W. the default

window is hamming

%window.

%06-Sep-04

%Esfandiar Zavarehei

if nargin<3

SP=.4;

end

if nargin<2

W=256;

end

if nargin<4

Window=hamming(W);

end

L=length(signal);

SP=fix(W.*SP);

N=fix((L-W)/SP +1); %number of segments

Index=(repmat(1:W,N,1)+repmat((0:(N-1))'*SP,1,W))';

hw=repmat(Window,1,N);

Seg=signal(Index).*hw;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [NoiseFlag, SpeechFlag, NoiseCounter,

Dist]=vad(signal,noise,NoiseCounter,NoiseMargin,Hangover)

%[NOISEFLAG, SPEECHFLAG, NOISECOUNTER,

DIST]=vad(SIGNAL,NOISE,NOISECOUNTER,NOISEMARGIN,HANGOVER)

%Spectral Distance Voice Activity Detector

%SIGNAL is the the current frames magnitude spectrum which is to be

labeld as %noise or speech, NOISE is noise magnitude spectrum template

(estimation), %NOISECOUNTER is the number of imediate previous noise

frames, NOISEMARGIN %(default 3)is the spectral distance threshold.

HANGOVER ( default 8 )is %the number of noise segments after which the

SPEECHFLAG is reset (goes to %zero). NOISEFLAG is set to one if the the

segment is labeld as noise %NOISECOUNTER returns the number of

previous noise segments, this value is %reset (to zero) whenever a speech

segment is detected. DIST is the %spectral distance.

%Note that:

% 1. This function works on a frame-by-frame basis not the whole % signal. If you would like to get the flags for the entire wave % signal you can use SEGMENT function and apply VAD to each % segment %2. It is safe to use the segment for updating noise when

%SPEECHFLAG is zero. NOISEFLAG would have been 1 for 8 frames %before that happens. %Author: Saeed Vaseghi

%edited by Esfandiar Zavarehei %Sep-04

if nargin<4

NoiseMargin=3;

end

if nargin<5

Hangover=8;

end

if nargin<3

NoiseCounter=0;

end

FreqResol=length(signal);

SpectralDist= 20*(log10(signal)-log10(noise));

SpectralDist(find(SpectralDist<0))=0;

Dist=mean(SpectralDist);

if (Dist < NoiseMargin)

NoiseFlag=1;

NoiseCounter=NoiseCounter+1; else

NoiseFlag=0;

NoiseCounter=0;

end

%Detect noise only periods and attenuate the signal if (NoiseCounter > Hangover) SpeechFlag=0;

else

SpeechFlag=1;

end

卡尔曼滤波算法与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函数参考 附录1.1 管理用命令 函数名功能描述函数名功能描述 addpath 增加一条搜索路径 rmpath 删除一条搜索路径 demo 运行Matlab演示程序 type 列出.M文件 doc 装入超文本文档 version 显示Matlab的版本号 help 启动联机帮助 what 列出当前目录下的有关文件 lasterr 显示最后一条信息 whatsnew 显示Matlab的新特性 lookfor 搜索关键词的帮助 which 造出函数与文件所在的目录 path 设置或查询Matlab路径 附录1.2管理变量与工作空间用命令 函数名功能描述函数名功能描述 clear 删除存中的变量与函数 pack 整理工作空间存 disp 显示矩阵与文本 save 将工作空间中的变量存盘 length 查询向量的维数 size 查询矩阵的维数 load 从文件中装入数据 who,whos 列出工作空间中的变量名 附录1.3文件与操作系统处理命令 函数名功能描述函数名功能描述 cd 改变当前工作目录 edit 编辑.M文件 delete 删除文件 matlabroot 获得Matlab的安装根目录 diary 将Matlab运行命令存盘 tempdir 获得系统的缓存目录 dir 列出当前目录的容 tempname 获得一个缓存(temp)文件 ! 执行操作系统命令 附录1.4窗口控制命令 函数名功能描述函数名功能描述 echo 显示文件中的Matlab中的命令 more 控制命令窗口的输出页面format 设置输出格式 附录1.5启动与退出命令 函数名功能描述函数名功能描述 matlabrc 启动主程序 quit 退出Matlab环境 startup Matlab自启动程序 附录2 运算符号与特殊字符附录 2.1运算符号与特殊字符 函数名功能描述函数名功能描述 + 加 ... 续行标志 - 减 , 分行符(该行结果不显示) * 矩阵乘 ; 分行符(该行结果显示) .* 向量乘 % 注释标志 ^ 矩阵乘方 ! 操作系统命令提示符 .^ 向量乘方矩阵转置 kron 矩阵kron积 . 向量转置 \ 矩阵左除 = 赋值运算 / 矩阵右除 == 关系运算之相等 .\ 向量左除 ~= 关系运算之不等 ./ 向量右除 < 关系运算之小于

matlab程序大全答案

频率特性类题目 1 一个系统的开环传递函数为 ，试绘制其当K=5、30时系统的开环频率特性Nyquist 图，并判断系统的稳定性。 w=linspace(0.5,5,1000)*pi; sys1=zpk([],[0 -10 -2],100) sys2=zpk([],[0 -10 -2],600) figure(1) nyquist(sys1,w) title('system nyquist charts with k=5') figure(2) nyquist(sys2,w) title('system nyquist charts with k=30') 由图可知K=5时，开环Nyquist 曲线没有包围（-1，j0）点，所以系统稳定。 K=30时，开环Nyquist 曲线包围（-1，j0）点，所以系统不稳定。 2系统开环传递函数为 ，建立其零极点增益模型， 然后分别绘制当K=5、K=30时系统的开环频率特性Bode 图，并判断系统的稳定性。 sys1=zpk([],[0 -10 -2],100) sys2=zpk([],[0 -10 -2],600) figure(1) [Gm1,Pm1,Wcg1,Wcp1]=margin(sys1) bode(sys1) title('system bode charts with k=5'),grid figure(2) [Gm2,Pm2,Wcg2,Wcp2]=margin(sys2) bode(sys2) title('system bode charts with k=30'),grid 因为K=5时，Wcg>Wcp,所以系统稳定。 K=10时，Wcg

扩展卡尔曼滤波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代码大全教学文案

m a t l a b代码大全

MATLAB主要命令汇总 MATLAB函数参考 附录1.1 管理用命令 函数名功能描述函数名功能描述 addpath 增加一条搜索路径 rmpath 删除一条搜索路径 demo 运行Matlab演示程序 type 列出.M文件 doc 装入超文本文档 version 显示Matlab的版本号 help 启动联机帮助 what 列出当前目录下的有关文件 lasterr 显示最后一条信息 whatsnew 显示Matlab的新特性 lookfor 搜索关键词的帮助 which 造出函数与文件所在的目录 path 设置或查询Matlab路径 附录1.2管理变量与工作空间用命令 函数名功能描述函数名功能描述 clear 删除内存中的变量与函数 pack 整理工作空间内存 disp 显示矩阵与文本 save 将工作空间中的变量存盘 length 查询向量的维数 size 查询矩阵的维数 load 从文件中装入数据 who,whos 列出工作空间中的变量名 附录1.3文件与操作系统处理命令 函数名功能描述函数名功能描述 cd 改变当前工作目录 edit 编辑.M文件 delete 删除文件 matlabroot 获得Matlab的安装根目录 diary 将Matlab运行命令存盘 tempdir 获得系统的缓存目录 dir 列出当前目录的内容 tempname 获得一个缓存(temp)文件 ! 执行操作系统命令 附录1.4窗口控制命令 函数名功能描述函数名功能描述 echo 显示文件中的Matlab中的命令 more 控制命令窗口的输出页面format 设置输出格式 附录1.5启动与退出命令 函数名功能描述函数名功能描述 matlabrc 启动主程序 quit 退出Matlab环境 startup Matlab自启动程序 附录2 运算符号与特殊字符附录

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

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

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

MATLAB简单程序大全

MATLAB简单程序大全 求特征值特征向量 A=[2 3 4;1 5 9;8 5 2] det(A) A' rank(A) inv(A) rref(A) eig(A)%求特征值和特征向量 卫星运行问题 h=200,H=51000,R=6378; a=(h+H+2*R)/2; c=(H-h)/2; b=(a^2-c^2)^(1/2); e=c/a; f=sqrt(1-exp(2).*cos(t)^2); l=int(f,t,0,pi/2) L=4*a.*l 动态玫瑰线 n=3;N=10000; theta=2*pi*(0:N)/N; r=cos(n*theta); x=r.*cos(theta); y=r.*sin(theta); comet(x,y) 二重积分 syms x y f=x^2*sin(y); int(int(f,x,0,1),y,0,pi) ezmesh(f,[0,1,0,pi]) 函数画图 syms x;f=exp(-0.2*x)*sin(0.5*x); ezplot(f,[0,8*pi])

玫瑰线 theta=0:0.01:2*pi; r=cos(3*theta); polar(theta,r,'r') 求x^2+y^2=1和x^2+z^2=1所围成的体积 syms x y z R r=1; Z=sqrt(1-x^2); y0=Z; V=8*int(int(Z,y,0,y0),x,0,1) 求导数及图像 f='1/(5+4*cos(x))'; subplot(1,2,1);ezplot(f) f1=diff(f) subplot(1,2,2);ezplot(f1) 绕x轴旋转 t=(0:20)*pi/10; r=exp(-.2*t).*sin(.5*t); theta=t; x=t'*ones(size(t)); y=r'*cos(theta); z=r'*sin(theta); mesh(x,y,z) colormap([0 0 0]) 某年是否闰年 year=input('input year:='); n1=year/4; n2=year/100; n3=year/400; if n1==fix(n1)&n2~=fix(n2) disp('是闰年') elseif n1==fix(n1)&n3==fix(n3) disp('是闰年') else

卡尔曼滤波入门简介及其算法MATLAB实现代码

卡尔曼滤波入门： 卡尔曼滤波是用来进行数据滤波用的，就是把含噪声的数据进行处理之后得出相对真值。卡尔曼滤波也可进行系统辨识。 卡尔曼滤波是一种基于统计学理论的算法，可以用来对含噪声数据进行在线处理，对噪声有特殊要求，也可以通过状态变量的增广形式实现系统辨识。 用上一个状态和当前状态的测量值来估计当前状态，这是因为上一个状态估计此时状态时会有误差，而测量的当前状态时也有一个测量误差，所以要根据这两个误差重新估计一个最接近真实状态的值。 信号处理的实际问题，常常是要解决在噪声中提取信号的问题，因此，我们需要寻找一种所谓有最佳线性过滤特性的滤波器。这种滤波器当信号与噪声同时输入时，在输出端能将信号尽可能精确地重现出来，而噪声却受到最大抑制。 维纳(Wiener)滤波与卡尔曼(Kalman)滤波就是用来解决这样一类从噪声中提取信号问题的一种过滤(或滤波)方法。 (1)过滤或滤波 - 从当前的和过去的观察值x(n)，x(n-1)，x(n-2)，…估计当前的信号值称为过滤或滤波; (2)预测或外推 - 从过去的观察值，估计当前的或将来的信号值称为预测或外推; (3)平滑或内插 - 从过去的观察值，估计过去的信号值称为平滑或内插; 因此，维纳过滤与卡尔曼过滤又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓“最佳”与“最优”是以最小均方误差为准则的。 维纳过滤与卡尔曼过滤都是解决最佳线性过滤和预测问题，并且都是以均方误差最小为准则的。因此在平稳条件下，它们所得到的稳态结果是一致的。然而，它们解决的方法有很大区别。 维纳过滤是根据全部过去的和当前的观察数据来估计信号的当前值，它的解是以均方误差最小条件下所得到的系统的传递函数H(z)或单位样本响应h(n)的形式给出的，因此更常称这种系统为最佳线性过滤器或滤波器。 而卡尔曼过滤是用前一个估计值和最近一个观察数据(它不需要全部过去的观察数据)来估计信号的当前值，它是用状态方程和递推的方法进行估计的，它的解是以估计值(常常是状态变量值)形式给出的。因此更常称这种系统为线性最优估计器或滤波器。 维纳滤波器只适用于平稳随机过程，而卡尔曼滤波器却没有这个限制。维纳过滤中信号和噪声是用相关函数表示的，因此设计维纳滤波器要求已知信号和噪声的相关函数。 卡尔曼过滤中信号和噪声是状态方程和量测方程表示的，因此设计卡尔曼滤波器要求已知状态方程和量测方程(当然，相关函数与状态方程和量测方程之间会存在一定的关系。卡尔曼过滤方法看来似乎比维纳过滤方法优越，它用递推法计算，不需要知道全部过去的数据，从而运用计算机计算方便，而且它可用于平稳和不平稳的随机过程(信号)，非时变和时变的系统。 但从发展历史上来看维纳过滤的思想是40年代初提出来的，1949年正式以书的形式出版。卡尔曼过滤到60年代初才提出来，它是在维纳过滤的基础上发展起来的，虽然如上所述它比维纳过滤方法有不少优越的地方，但是最佳线性过滤问题是由维纳过滤首先解决的，维纳过滤的物理概念比较清楚，也可以认为卡尔曼滤波仅仅是对最佳线性过滤问题提出的一种新的算法。 卡尔曼滤波在数学上是一种统计估算方法，通过处理一系列带有误差的实际量测数据而得到的物理参数的最佳估算。例如在气象应用上，根据滤波的基本思想，利用前一时刻预报误差的反馈信息及时修正预报方程，以提高下一时刻预报精度。作温度预报一般只需要连续两个月的资料即可建立方程和递推关系。

matlab的编码大全

附录Matlab源程序 附录A 信息熵 % 函数说明：% % H=entropy(P,r) 为信息熵函数% % P为信源的概率矢量, r为进制数% % H为信息熵% %****************************** % function H=entropy(P,r) if (length(find(P<=0))~=0) error('Not a prob.vector,negative component'); % 判断是否符合概率分布条件end if (abs(sum(P)-1)>10e-10) error('Not a prob.vector,component do not add up to 1'); end H=(sum(-P.*log2(P)))/(log2(r)+eps); 附录B 离散无记忆信道容量的迭代计算 % 信道容量C的迭代算法% % 函数说明：% % [CC,Paa]=ChannelCap(P,k) 为信道容量函数% % 变量说明：% % P：输入的正向转移概率矩阵，k：迭代计算精度% % CC：最佳信道容量，Paa：最佳输入概率矩阵% % Pa：初始输入概率矩阵，Pba：正向转移概率矩阵% % Pb：输出概率矩阵，Pab：反向转移概率矩阵% % C：初始信道容量，r：输入符号数，s：输出符号数% %************************************************** % function [CC,Paa]=ChannelCap(P,k) % 提示错误信息 if (length(find(P<0)) ~=0) error('Not a prob.vector,negative component'); % 判断是否符合概率分布条件end

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

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

matlab快捷键大全

1.使用“Tab”键可以在Matlab进行命令输入时补全变量名或者命令名； 2.使用“Ctrl+C”组合键能够强制从运行的或者进入死循环的Matlab程序中退出； 3.使用“上”“下”方向键能够调用Matlab的历史命令； 4.Matlab的变量命名规则注意： a.区分大小写； b.必须以字母开头； c.中间不能有空格、标点符号等； 5.在一条命令或语句的末尾加上一个分号，则Matlab就不会在屏幕上显示这条命令计算的结果； 6.在Matlab中，“...”(3个句点)称为“续行符”，表示同一语句的延续输入； 注意：只有当续行符出现在变量名和运算符之间，才能起到语句延续的作用。出现在一个变量名中间时，是不能实现语句延续的。换句话说，变量名不能分割成两行书写。 7.在进行三角函数运算时，Matlab使用的是“弧度制”而不是“角度制”； 8.选择需要注释的语句块以后，使用“Ctrl+R”可以进行Matlab语句注释；使用“Ctrl+T”可以进行语句块的解注释； 9.Matlab数组的下标是从“1”开始的，这点与C语言不同； 1 0.使用空格和逗号隔开的元素指定的是同一行的元素，使用分号或者回车分开的元素指定的是不同行的元素。 Q1：matlab有没有监视内存的方法? A：用函数whos。 Q2：如何解决matlab7.0命令窗口跳出一大堆java错误...

A：换matlab7的sp2。 Q3：自从安装matlab，一开机就在进程里有matlab。 能不能开机的时候进程就不运行matlab？ A：开始-->控制面板-->管理工具-->服务 把MATLAB Server的属性改成“手动”就行了。 Q4：退出matlab7程序运行的快捷键。 A：ctrl+q Q5：matlab7远程控制是否有限制？ A：不能远程控制,不过你可以先在你的remote机器上打开,然后就可以用了。 Q6：Matlab占用资源太多怎么办？ A：用matlab-nojvm启动（如果不需要图形界面）。 Q7：怎样给matlab添加新的toolbox？ A：在matlab的文件菜单里边添加路径，选set path。 Q8：请问matlab有没有命令可以调出历史输入啊？ A：command history窗口。 Q9：matlab7.0不能在64位的cpu下运行？ A：matlab应该是依赖于自己的虚拟机的 但是好像这个虚拟机是在IA32里面作出来的，所以，应该找个带64位的java虚拟机替换原来的，不过不一定能行or so,记不清了） Q10：matlab有没有注释一段的功能？ A：选中一段代码，ctrl r就是区段注释 选中一段代码，ctrl t取消区段注释

matlab 指令大全

分享 我的分享 当前分享 返回分享首页? 分享 matlab命令，应该很全了！来源：李家叶的日志 matlab命令 一、常用对象操作：除了一般windows窗口的常用功能键外。 1、!dir 可以查看当前工作目录的文件。!dir& 可以在dos状态下查看。 2、who 可以查看当前工作空间变量名，whos 可以查看变量名细节。 3、功能键： 功能键快捷键说明 方向上键Ctrl+P 返回前一行输入 方向下键Ctrl+N 返回下一行输入 方向左键Ctrl+B 光标向后移一个字符 方向右键Ctrl+F 光标向前移一个字符 Ctrl+方向右键Ctrl+R 光标向右移一个字符 Ctrl+方向左键Ctrl+L 光标向左移一个字符 home Ctrl+A 光标移到行首 End Ctrl+E 光标移到行尾 Esc Ctrl+U 清除一行 Del Ctrl+D 清除光标所在的字符 Backspace Ctrl+H 删除光标前一个字符 Ctrl+K 删除到行尾 Ctrl+C 中断正在执行的命令 4、clc可以命令窗口显示的内容，但并不清除工作空间。 二、函数及运算 1、运算符： ＋：加，－：减，*：乘，/：除，\：左除^：幂，‘：复数的共轭转置，（）：制定运算顺序。 2、常用函数表： sin( ) 正弦（变量为弧度） Cot( ) 余切（变量为弧度） sind( ) 正弦（变量为度数） Cotd( ) 余切（变量为度数） asin( ) 反正弦（返回弧度） acot( ) 反余切（返回弧度） Asind( ) 反正弦（返回度数） acotd( ) 反余切（返回度数） cos( ) 余弦（变量为弧度） exp( ) 指数 cosd( ) 余弦（变量为度数）

(完整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;

实验一 典型环节的MATLAB仿真汇总

实验一 典型环节的MATLAB 仿真 一、实验目的 1．熟悉MATLAB 桌面和命令窗口，初步了解SIMULINK 功能模块的使用方法。 2．通过观察典型环节在单位阶跃信号作用下的动态特性，加深对各典型环节响应曲线的理解。 3．定性了解各参数变化对典型环节动态特性的影响。 二、SIMULINK 的使用 MATLAB 中SIMULINK 是一个用来对动态系统进行建模、仿真和分析的软件包。利用SIMULINK 功能模块可以快速的建立控制系统的模型，进行仿真和调试。 1．运行MATLAB 软件，在命令窗口栏“>>”提示符下键入simulink 命令，按Enter 键或在工具栏单击按钮，即可进入如图1-1所示的SIMULINK 仿真 环境下。 2．选择File 菜单下New 下的Model 命令，新建一个simulink 仿真环境常规模板。 3．在simulink 仿真环境下，创建所需要的系统 三、实验内容 按下列各典型环节的传递函数，建立相应的SIMULINK 仿真模型，观察并记录其单位阶跃响应波形。 ① 比例环节1)(1=s G 和2)(1=s G 实验处理：1)(1=s G SIMULINK 仿真模型

波形图为： 实验处理：2)(1=s G SIMULINK 仿真模型 波形图为： 实验结果分析：增加比例函数环节以后，系统的输出型号将输入信号成倍数放大. ② 惯性环节11)(1+= s s G 和15.01)(2+=s s G 实验处理：1 1 )(1+=s s G SIMULINK 仿真模型

波形图为： 实验处理：1 5.01 )(2+= s s G SIMULINK 仿真模型 波形图为： 实验结果分析：当1 1 )(1+= s s G 时，系统达到稳定需要时间接近5s,当

卡尔曼滤波简介及其算法实现代码

卡尔曼滤波简介及其算法实现代码 卡尔曼滤波算法实现代码（C，C＋＋分别实现） 卡尔曼滤波器简介 近来发现有些问题很多人都很感兴趣。所以在这里希望能尽自己能力跟大家讨论一些力所能及的算法。现在先讨论一下卡尔曼滤波器，如果时间和能力允许，我还希望能够写写其他的算法，例如遗传算法，傅立叶变换，数字滤波，神经网络，图像处理等等。 因为这里不能写复杂的数学公式，所以也只能形象的描述。希望如果哪位是这方面的专家，欢迎讨论更正。 卡尔曼滤波器– Kalman Filter 1．什么是卡尔曼滤波器 （What is the Kalman Filter?） 在学习卡尔曼滤波器之前，首先看看为什么叫“卡尔曼”。跟其他著名的理论（例如傅立叶变换，泰勒级数等等）一样，卡尔曼也是一个人的名字，而跟他们不同的是，他是个现代人！ 卡尔曼全名Rudolf Emil Kalman，匈牙利数学家，1930年出生于匈牙利首都布达佩斯。1953，1954年于麻省理工学院分别获得电机工程学士及硕士学位。1957年于哥伦比亚大学获得博士学位。我们现在要学习的卡尔曼滤波器，正是源于他的博士论文和1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》（线性滤波与预测问题的新方法）。如果对这编论文有兴趣，可以到这里的地址下载： https://www.doczj.com/doc/f210740018.html,/~welch/media/pdf/Kalman1960.pdf。 简单来说，卡尔曼滤波器是一个“optimal recursive data processing algorithm（最优化自回归数据处理算法）”。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的。他的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 2．卡尔曼滤波器的介绍 （Introduction to the Kalman Filter） 为了可以更加容易的理解卡尔曼滤波器，这里会应用形象的描述方法来讲解，而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是，他的5条公式是其核心内容。结合现代的计算机，其实卡尔曼的程序相当的简单，只要你理解了他的那5条公式。 在介绍他的5条公式之前，先让我们来根据下面的例子一步一步的探索。 假设我们要研究的对象是一个房间的温度。根据你的经验判断，这个房间的温度是恒定的，也就

数字图像处理代码大全

1.图像反转 MATLAB程序实现如下： I=imread('xian.bmp'); J=double(I); J=-J+(256-1); %图像反转线性变换 H=uint8(J); subplot(1,2,1),imshow(I); subplot(1,2,2),imshow(H); 2.灰度线性变换 MATLAB程序实现如下： I=imread('xian.bmp'); subplot(2,2,1),imshow(I); title('原始图像'); axis([50,250,50,200]); axis on; %显示坐标系 I1=rgb2gray(I); subplot(2,2,2),imshow(I1); title('灰度图像'); axis([50,250,50,200]); axis on; %显示坐标系 J=imadjust(I1,[0.1 0.5],[]); %局部拉伸，把[0.1 0.5]的灰度拉伸为[0 1]

subplot(2,2,3),imshow(J); title('线性变换图像[0.1 0.5]'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 K=imadjust(I1,[0.3 0.7],[]); %局部拉伸，把[0.3 0.7]的灰度拉伸为[0 1] subplot(2,2,4),imshow(K); title('线性变换图像[0.3 0.7]'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 3.非线性变换 MATLAB程序实现如下： I=imread('xian.bmp'); I1=rgb2gray(I); subplot(1,2,1),imshow(I1); title('灰度图像'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 J=double(I1);

卡尔曼滤波两例题含matlab程序汇总

设高度的测量误差是均值为0、方差为1的高斯白噪声随机序列，该物体的初始高度0h 和速度0V 也是高斯分布的随机变量，且0000019001000,var 10/02Eh h m P EV m s V ???????? ===? ??????? ???? ????。试求该物体高度和速度随时间变化的最优估计。（2/80.9s m g =） 解： 1. 令()()()h k X k v k ?? =? ??? t=1 R （k ）=1 Q(k)=0 根据离散时间卡尔曼滤波公式，则有： (1)(1,)()()X k k k X k U k φ+=++ (1)(1)(1)(1)Y k H k X k V k +=++++ (1,)k k φ+＝ 11t -?? ??? ? ()U k ＝ 20.5gt gt ??-???? (1)H k +＝[]10 滤波初值：^ 1900(0|0)(0)10X EX ?? ==???? 0100(0|0)var[(0)]2P X P ?? ===? ??? 一步预测：^^ (1|)(1,)(|)()X k k k k X k k U k φ+=++ (1|)(1,)(|)(1,)T P k k k k P k k k k φφ+=++ 滤波增益：1 (1)(1|)(1)[(1)(1|)(1)(1)]T T K k P k k H k H k P k k H k R k -+=+++++++ 滤波计算：^ ^ ^ (1|1)(1|)(1)[(1)(1)(1|)]X k k X k k K k Y k H k X k k ++=++++-++ (1|1)[(1)(1)](1|)P k k I K k H k P k k ++=-+++ 2. 实验结果

Matlab代码1运行结果

when rho=0.7000 epsilon_c2=0.9760183979. when rho=0.6500 epsilon_c1=0.9545430725. when rho=0.6500 epsilon_c2=0.9808187253. when rho=0.6000 epsilon_c1=0.9637291916. when rho=0.6000 epsilon_c2=0.9844500729. when rho=0.5500 epsilon_c1=0.9709578854. when rho=0.5500 epsilon_c2=0.9872930980. when rho=0.5000 epsilon_c1=0.9767709652. when rho=0.5000 epsilon_c2=0.9895793686. when rho=0.4500 epsilon_c1=0.9815214506. when rho=0.4500 epsilon_c2=0.9914578443. when rho=0.4000 epsilon_c1=0.9854483433. when rho=0.4000 epsilon_c2=0.9930286972. when rho=0.3500 epsilon_c1=0.9887183111. when rho=0.3500 epsilon_c2=0.9943617736. when rho=0.3000 epsilon_c1=0.9914501549. when rho=0.3000 epsilon_c2=0.9955072624. when rho=0.2500 epsilon_c1=0.9937297787. when rho=0.2500 epsilon_c2=0.9965021564. when rho=0.2000 epsilon_c1=0.9956196526. when rho=0.2000 epsilon_c2=0.9973743212. when rho=0.1500 epsilon_c1=0.9971649360. when rho=0.1500 epsilon_c2=0.9981451473. when rho=0.1000 epsilon_c1=0.9983974871. when rho=0.1000 epsilon_c2=0.9988313269. when rho=0.0500 epsilon_c1=0.9993384730. when rho=0.0500 epsilon_c2=0.9994460773. when rho=0.0000 epsilon_c1=1.0000000000. when rho=0.0000 epsilon_c2=1.0000000000. out = Columns 1 through 11 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.0000 0.6540 0.8011 0.8661 0.9027 0.9262 0.9425 0.9545 0.9637 0.9710 0.9768 0 0.8333 0.9134 0.9437 0.9596 0.9694 0.9760 0.9808 0.9845 0.9873 0.9896 Columns 12 through 21 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000