搜档网
当前位置:搜档网 › KSVD原子库训练MATLAB程序

KSVD原子库训练MATLAB程序

KSVD原子库训练MATLAB程序
KSVD原子库训练MATLAB程序

% KSVD running file

% in this file a synthetic test of the K-SVD algorithm is performed. First,

% a random dictionary with normalized columns is being generated, and then

% a set of data signals, each as a linear combination of 3 dictionary

% element is created, with noise level of 20SNR. this set is given as input

% to the K-SVD algorithm.

% a different mode for activating the K-SVD algorithm is until a fixed

% error is reached in the Sparse coding stage, instead until a fixed number of coefficients is found

% (it was used by us for the

% denoising experiments). in order to switch between those two modes just

% change the param.errorFlag (0 - for fixed number of coefficients, 1 -

% until a certain error is reached).

param.L = 3; % number of elements in each linear combination.

param.K = 50; % number of dictionary elements

param.numIteration = 50; % number of iteration to execute the K-SVD algorithm. param.errorFlag = 0; % decompose signals until a certain error is reached. do not use fix number of coefficients.

%param.errorGoal = sigma;

param.preserveDCAtom = 0;

%%%%%%% creating the data to train on %%%%%%%%

N = 1500; % number of signals to generate

n = 20; % dimension of each data

SNRdB = 20; % level of noise to be added

[param.TrueDictionary, D, x] = gererateSyntheticDictionaryAndData(N, param.L, n, param.K, SNRdB); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%

%%%%%%%% initial dictionary: Dictionary elements %%%%%%%%

param.InitializationMethod = 'DataElements';

param.displayProgress = 1;

disp('Starting to train the dictionary');

[Dictionary,output] = KSVD(D,param);

disp(['The KSVD algorithm retrived ',num2str(output.ratio(end)),' atoms from the original dictionary']);

[Dictionary,output] = MOD(D,param);

disp(['The MOD algorithm retrived ',num2str(output.ratio(end)),' atoms from the original dictionary']);

function [Dictionary, data, coefs] = gererateSyntheticDictionaryAndData(N, L, dim, K, SNRdB)

randn('state',sum(100*clock));

rand('state',sum(100*clock));

Dictionary = randn(dim,K);

Dictionary = Dictionary*diag(1./sqrt(sum(Dictionary.*Dictionary)));

[data,coefs] = CreateDataFromDictionarySimple(Dictionary, N, L);

if (SNRdB==0) | (SNRdB == 80)

return

else

noise = randn(size(data));

actualNoise = calcNoiseFromSNR(SNRdB,data, noise);

SNR = calcSNR(data, data+actualNoise);

data = data + actualNoise*SNR/SNRdB;

end

function [D,xOrig] = CreateDataFromDictionarySimple(dictionary, numElements, numCoef)

maxRangeOfCoef = 1;

resolution = 0.0001;

xOrig = zeros(size(dictionary,2),numElements);

%vecOfValues = -1*maxRangeOfCoef:resolution:maxRangeOfCoef;

%coefs = randsrc(numCoef,numElements,vecOfValues);

coefs = randn(numCoef,numElements)*maxRangeOfCoef;

xOrig(1:numCoef,:) = coefs;

for i=1:size(xOrig,2)

xOrig(:,i) = xOrig(randperm(size(xOrig,1)),i);

end

%dictionaryElementIndices =

randsrc(numCoef*numElements,1,[1:size(dictionary,2)]) ;

%matrixOfIndices = repmat([1:numElements],numCoef,1);

%xOrig(sub2ind(size(xOrig),dictionaryElementIndices,matrixOfIndices(:))) = coefs;

D = dictionary*xOrig;

function actualNoise = calcNoiseFromSNR(TargerSNR, signal, randomNoise)

signal = signal(:);

randomNoiseRow = randomNoise(:);

signal_2 = sum(signal.^2);

ActualNoise_2 = signal_2/(10^(TargerSNR/10));

noise_2 = sum(randomNoiseRow.^2);

ratio = ActualNoise_2./noise_2;

actualNoise = randomNoiseRow.*repmat(sqrt(ratio),size(randomNoiseRow,1),1); actualNoise = reshape(actualNoise,size(randomNoise));

function SNR = calcSNR(origSignal, noisySignal)

errorSignal = origSignal-noisySignal;

signal_2 = sum(origSignal.^2);

noise_2 = sum(errorSignal.^2);

SNRValues = 10*log10(signal_2./noise_2);

SNR = mean(SNRValues);

function [Dictionary,output] = KSVD(...

Data,... % an nXN matrix that contins N signals (Y), each of dimension n.

param)

%

=============================================================== ==========

% K-SVD algorithm

%

=============================================================== ==========

% The K-SVD algorithm finds a dictionary for linear representation of

% signals. Given a set of signals, it searches for the best dictionary that

% can sparsely represent each signal. Detailed discussion on the algorithm

% and possible applications can be found in "The K-SVD: An Algorithm for

% Designing of Overcomplete Dictionaries for Sparse Representation", written

% by M. Aharon, M. Elad, and A.M. Bruckstein and appeared in the IEEE Trans.

% On Signal Processing, Vol. 54, no. 11, pp. 4311-4322, November 2006.

%

=============================================================== ==========

% INPUT ARGUMENTS:

% Data an nXN matrix that contins N signals (Y), each of dimension n.

% param structure that includes all required

% parameters for the K-SVD execution.

% Required fields are:

% K, ... the number of dictionary elements to train

% numIteration,... number of iterations to perform.

% errorFlag... if =0, a fix number of coefficients is

% used for representation of each signal. If so, param.L must be

% specified as the number of representing atom. if =1, arbitrary number

% of atoms represent each signal, until a specific representation error

% is reached. If so, param.errorGoal must be specified as the allowed

% error.

% preserveDCAtom... if =1 then the first atom in the dictionary

% is set to be constant, and does not ever change. This

% might be useful for working with natural

% images (in this case, only param.K-1

% atoms are trained).

% (optional, see errorFlag) L,... % maximum coefficients to use in OMP coefficient calculations.

% (optional, see errorFlag) errorGoal, ... % allowed representation error in representing each signal.

% InitializationMethod,... mehtod to initialize the dictionary, can

% be one of the following arguments:

% * 'DataElements' (initialization by the signals themselves), or:

% * 'GivenMatrix' (initialization by a given matrix param.initialDictionary).

% (optional, see InitializationMethod) initialDictionary,... % if the initialization method

% is 'GivenMatrix', this is the matrix that will be used.

% (optional) TrueDictionary, ... % if specified, in each

% iteration the difference between this dictionary and the trained one

% is measured and displayed.

% displayProgress, ... if =1 progress information is displyed. If

param.errorFlag==0,

% the average repersentation error (RMSE) is displayed, while if

% param.errorFlag==1, the average number of required coefficients for

% representation of each signal is displayed.

%

=============================================================== ==========

% OUTPUT ARGUMENTS:

% Dictionary The extracted dictionary of size nX(param.K).

% output Struct that contains information about the current

run. It may include the following fields:

% CoefMatrix The final coefficients matrix (it should hold that Data equals approximately Dictionary*output.CoefMatrix.

% ratio If the true dictionary was defined (in

% synthetic experiments), this parameter holds a vector of length

% param.numIteration that includes the detection ratios in each

% iteration).

% totalerr The total representation error after each

% iteration (defined only if

% param.displayProgress=1 and

% param.errorFlag = 0)

% numCoef A vector of length param.numIteration that

% include the average number of coefficients required for representation

% of each signal (in each iteration) (defined only if

% param.displayProgress=1 and

% param.errorFlag = 1)

%

=============================================================== ==========

if (~isfield(param,'displayProgress'))

param.displayProgress = 0;

end

totalerr(1) = 99999;

if (isfield(param,'errorFlag')==0)

param.errorFlag = 0;

end

if (isfield(param,'TrueDictionary'))

displayErrorWithTrueDictionary = 1;

ErrorBetweenDictionaries = zeros(param.numIteration+1,1);

ratio = zeros(param.numIteration+1,1);

else

displayErrorWithTrueDictionary = 0;

ratio = 0;

end

if (param.preserveDCAtom>0)

FixedDictionaryElement(1:size(Data,1),1) = 1/sqrt(size(Data,1));

else

FixedDictionaryElement = [];

end

% coefficient calculation method is OMP with fixed number of coefficients

if (size(Data,2) < param.K)

disp('Size of data is smaller than the dictionary size. Trivial solution...');

Dictionary = Data(:,1:size(Data,2));

return;

elseif (strcmp(param.InitializationMethod,'DataElements'))

Dictionary(:,1:param.K-param.preserveDCAtom) =

Data(:,1:param.K-param.preserveDCAtom);

elseif (strcmp(param.InitializationMethod,'GivenMatrix'))

Dictionary(:,1:param.K-param.preserveDCAtom) =

param.initialDictionary(:,1:param.K-param.preserveDCAtom);

end

% reduce the components in Dictionary that are spanned by the fixed

% elements

if (param.preserveDCAtom)

tmpMat = FixedDictionaryElement \ Dictionary;

Dictionary = Dictionary - FixedDictionaryElement*tmpMat;

end

%normalize the dictionary.

Dictionary = Dictionary*diag(1./sqrt(sum(Dictionary.*Dictionary)));

Dictionary = Dictionary.*repmat(sign(Dictionary(1,:)),size(Dictionary,1),1); % multiply in the sign of the first element.

totalErr = zeros(1,param.numIteration);

% the K-SVD algorithm starts here.

for iterNum = 1:param.numIteration

% find the coefficients

if (param.errorFlag==0)

%CoefMatrix = mexOMPIterative2(Data, [FixedDictionaryElement,Dictionary],param.L);

CoefMatrix = OMP([FixedDictionaryElement,Dictionary],Data, param.L);

else

%CoefMatrix = mexOMPerrIterative(Data, [FixedDictionaryElement,Dictionary],param.errorGoal);

CoefMatrix = OMPerr([FixedDictionaryElement,Dictionary],Data,

param.errorGoal);

param.L = 1;

end

replacedVectorCounter = 0;

rPerm = randperm(size(Dictionary,2));

for j = rPerm

[betterDictionaryElement,CoefMatrix,addedNewVector] =

I_findBetterDictionaryElement(Data,...

[FixedDictionaryElement,Dictionary],j+size(FixedDictionaryElement,2),...

CoefMatrix ,param.L);

Dictionary(:,j) = betterDictionaryElement;

if (param.preserveDCAtom)

tmpCoef = FixedDictionaryElement\betterDictionaryElement;

Dictionary(:,j) = betterDictionaryElement -

FixedDictionaryElement*tmpCoef;

Dictionary(:,j) = Dictionary(:,j)./sqrt(Dictionary(:,j)'*Dictionary(:,j));

end

replacedVectorCounter = replacedVectorCounter+addedNewVector;

end

if (iterNum>1 & param.displayProgress)

if (param.errorFlag==0)

output.totalerr(iterNum-1) =

sqrt(sum(sum((Data-[FixedDictionaryElement,Dictionary]*CoefMatrix).^2))/prod(size(D ata)));

disp(['Iteration ',num2str(iterNum),' Total error is:

',num2str(output.totalerr(iterNum-1))]);

else

output.numCoef(iterNum-1) = length(find(CoefMatrix))/size(Data,2);

disp(['Iteration ',num2str(iterNum),' Average number of coefficients: ',num2str(output.numCoef(iterNum-1))]);

end

end

if (displayErrorWithTrueDictionary )

[ratio(iterNum+1),ErrorBetweenDictionaries(iterNum+1)] =

I_findDistanseBetweenDictionaries(param.TrueDictionary,Dictionary);

disp(strcat(['Iteration ', num2str(iterNum),' ratio of restored elements:

',num2str(ratio(iterNum+1))]));

output.ratio = ratio;

end

Dictionary =

I_clearDictionary(Dictionary,CoefMatrix(size(FixedDictionaryElement,2)+1:end,:),Data) ;

if (isfield(param,'waitBarHandle'))

waitbar(iterNum/param.counterForWaitBar);

end

end

output.CoefMatrix = CoefMatrix;

Dictionary = [FixedDictionaryElement,Dictionary]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% findBetterDictionaryElement %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [betterDictionaryElement,CoefMatrix,NewVectorAdded] =

I_findBetterDictionaryElement(Data,Dictionary,j,CoefMatrix,numCoefUsed)

if (length(who('numCoefUsed'))==0)

numCoefUsed = 1;

end

relevantDataIndices = find(CoefMatrix(j,:)); % the data indices that uses the j'th dictionary element.

if (length(relevantDataIndices)<1) %(length(relevantDataIndices)==0) ErrorMat = Data-Dictionary*CoefMatrix;

ErrorNormVec = sum(ErrorMat.^2);

[d,i] = max(ErrorNormVec);

betterDictionaryElement = Data(:,i);%ErrorMat(:,i); %

betterDictionaryElement =

betterDictionaryElement./sqrt(betterDictionaryElement'*betterDictionaryElement);

betterDictionaryElement =

betterDictionaryElement.*sign(betterDictionaryElement(1));

CoefMatrix(j,:) = 0;

NewVectorAdded = 1;

return;

end

NewVectorAdded = 0;

tmpCoefMatrix = CoefMatrix(:,relevantDataIndices);

tmpCoefMatrix(j,:) = 0;% the coeffitients of the element we now improve are not relevant.

errors =(Data(:,relevantDataIndices) - Dictionary*tmpCoefMatrix); % vector of errors that we want to minimize with the new element

% % the better dictionary element and the values of beta are found using svd.

% % This is because we would like to minimize || errors - beta*element ||_F^2.

% % that is, to approximate the matrix 'errors' with a one-rank matrix. This

% % is done using the largest singular value.

[betterDictionaryElement,singularValue,betaVector] = svds(errors,1);

CoefMatrix(j,relevantDataIndices) = singularValue*betaVector';% *signOfFirstElem %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% findDistanseBetweenDictionaries %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ratio,totalDistances] = I_findDistanseBetweenDictionaries(original,new)

% first, all the column in oiginal starts with positive values.

catchCounter = 0;

totalDistances = 0;

for i = 1:size(new,2)

new(:,i) = sign(new(1,i))*new(:,i);

end

for i = 1:size(original,2)

d = sign(original(1,i))*original(:,i);

distances =sum ( (new-repmat(d,1,size(new,2))).^2);

[minValue,index] = min(distances);

errorOfElement = 1-abs(new(:,index)'*d);

totalDistances = totalDistances+errorOfElement;

catchCounter = catchCounter+(errorOfElement<0.01);

end

ratio = 100*catchCounter/size(original,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % I_clearDictionary %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Dictionary = I_clearDictionary(Dictionary,CoefMatrix,Data)

T2 = 0.99;

T1 = 3;

K=size(Dictionary,2);

Er=sum((Data-Dictionary*CoefMatrix).^2,1); % remove identical atoms

G=Dictionary'*Dictionary; G = G-diag(diag(G));

for jj=1:1:K,

if max(G(jj,:))>T2 | length(find(abs(CoefMatrix(jj,:))>1e-7))<=T1 , [val,pos]=max(Er);

Er(pos(1))=0;

Dictionary(:,jj)=Data(:,pos(1))/norm(Data(:,pos(1)));

G=Dictionary'*Dictionary; G = G-diag(diag(G));

end;

end;

MATLAB编程作业

《Matlab 编程训练》 作业 专 业 学生姓名 班级 学 号 指导教师 完成日期

实训一 MATLAB 语言介绍和数值计算 1.先求下列表达式的值,然后显示MATLAB 工作空间的使用情况并保存变量。 12 2sin851z e =+ . 2. 已知 1234413134787,2033657327A B --???? ????==???? ????-???? ,求下列表达式的值: (1) A+6*B 和A-B+I (其中I 为单位矩阵) A+6*B:

A-B+I: (2)A*B和A.*B A*B程序: A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] c=A*B 结果: A.*B程序: A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] D=A.*B 结果:

(3)A^3和A.^3 A^3程序: A=[12 34 -4;34 7 87;3 65 7] E=A^3 结果: A.^3程序: A=[12 34 -4;34 7 87;3 65 7] C=A.^3 (4)A/B及B\A A/B程序: A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] C=A/B 结果:

B\A程序: A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] D=B\A 结果: (5)将矩阵C=B\A的右下角2*2子矩阵赋给D, 并(3)保存变量(mat文件)程序: A=[12 34 -4;34 7 87;3 65 7]; B=[1 3 -1;2 0 3;3 -2 7]; C=B*inv(A); D=C(2:3,2:3) 结果:

MATLAB程序设计教程(第二版)课后答案(可编辑修改word版)

MATLAB 第二版课后答案unit3-8 unit3 实验指导 1、 n=input('请输入一个三位数:'); a=fix(n/100); b=fix((n-a*100)/10); c=n-a*100-b*10; d=c*100+b*10+a 2(1) n=input('请输入成绩'); switch n case num2cell(90:100) p='A'; case num2cell(80:89) p='B'; case num2cell(70:79) p='C'; case num2cell(60:69) p='D'; otherwise p='E'; end price=p (2)n=input('请输入成绩'); if n>=90&n<=100 p='A'; elseif n>=80&n<=89 p='B'; elseif n>=70&n<=79 p='C'; elseif n>=60&n<=69 p='D'; else p='E'; end price=p (3)try n; catch price='erroe' end 3 n=[1,5,56,4,3,476,45,6,3,76,45,6,4,3,6,4,23,76,908,6];

b=n(1); for m=2:20 if n(m)>a a=n(m); elseif n(m)=0 disp(A([n],:)); elseif n<0 disp(lasterr); else disp(A([6],:)); disp(lasterr); end 7(1) f=[];

实验二 Matlab程序设计基本方法1

实验二Matlab程序设计基本方法 覃照乘自092 电气工程学院 一、实验目的: 1、熟悉MATLAB 程序编辑与设计环境 2、掌握各种编程语句语法规则及程序设计方法 3、函数文件的编写和设计 4、了解和熟悉跨空间变量传递和赋值 二、实验基本知识: ◆for循环结构 语法:for i=初值:增量:终值 语句1 …… 语句n end 说明:1.i=初值:终值,则增量为1。 2.初值、增量、终值可正可负,可以是整数,也可以是小数,只须符合数学逻辑。 ◆while 循环结构 语法:while 逻辑表达式 循环体语句 end 说明:1、whiIe结构依据逻辑表达式的值判断是否执行循环体语勾。若表达式的值为真,执行循环体语句一次、在反复执行时,每次都要进行判断。若表达 式的值为假,则程序执行end之后的语句。 2、为了避免因逻辑上的失误,而陷入死循环,建议在循环体语句的适当位 置加break语句、以便程序能正常执行。(执行循环体的次数不确定; 每一次执行循环体后,一定会改变while后面所跟关系式的值。) 3、while循环也可以嵌套、其结构如下:

while逻辑表达式1 循环体语句1 while逻辑表达式2 循环体语句2 end 循环体语句3 end ◆if-else-end分支结构 if 表达式1 语句1 else if 表达式2(可选) 语句2 else(可选) 语句3 end end 说明:1.if结构是一个条件分支语句,若满足表达式的条件,则往下执行;若不满足,则跳出if结构。 2.else if表达式2与else为可选项,这两条语句可依据具体情况取舍。 3.注意:每一个if都对应一个end,即有几个if,记就应有几个end。 ◆switch-case结构 语法:switch表达式 case常量表达式1 语句组1 case常量表达式2 语句组2 …… otherwise 语句组n end

MATLAB程序设计作业

Matlab程序设计 班级 姓名 学号

《MATLAB程序设计》作业 1、考虑如下x-y 一组实验数据: x=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10] y=[1.2, 3, 4, 4, 5, 4.7, 5, 5.2, 6, 7.2] 分别绘出plot的原始数据、一次拟合曲线和三次拟合曲线,给出MATLAB代码和运行结果。 代码如下: x=[1,2,3,4,5,6,7,8,9,10]; y=[1.2,3,4,4,5,4.7,5,5.2,6,7.2]; plot(x,y); title('原始数据'); p=polyfit(x,y,1); q=polyval(p,x); figure,plot(x,q); title('一次拟合'); p=polyfit(x,y,2); q=polyval(p,x); figure,plot(x,q); title('二次拟合'); 运行结果如下:

1 2 3 4 5 6 7 8 9 10 12 3 4 5 6 7 8 原始数据 123 456789 102 2.5 3 3.54 4.5 55.56 6.57一次拟合 123456789 101 2 3 4 5 6 7 二次拟合 2、在[0,3π]区间,绘制y=sin(x)曲线(要求消去负半波,即(π,2π)区间内的函数值置零),求出曲线y 的平均值,以及y 的最大值及其最大值的位置。给出执行代码和运行结果。 代码如下: clear clc x=(0:0.01:3*pi); y=sin(x); plot(x,y); y1=(y>=0).*y; figure,plot(x,y1);

MATLAB程序设计教程课后答案

实验指导 1、 n=input('请输入一个三位数:'); a=fix(n/100); b=fix((n-a*100)/10); c=n-a*100-b*10; d=c*100+b*10+a 2(1) n=input('请输入成绩'); switch n case num2cell(90:100) p='A'; case num2cell(80:89) p='B'; case num2cell(70:79) p='C'; case num2cell(60:69) p='D'; otherwise p='E'; end price=p (2)n=input('请输入成绩'); if n>=90&n<=100 p='A'; elseif n>=80&n<=89 p='B'; elseif n>=70&n<=79 p='C'; elseif n>=60&n<=69 p='D'; else p='E'; end price=p (3)try n; catch price='erroe' end 3 n=[1,5,56,4,3,476,45,6,3,76,45,6,4,3,6,4,23,76,908,6]; a=n(1);

for m=2:20 if n(m)>a a=n(m); elseif n(m)=0 disp(A([n],:)); elseif n<0 disp(lasterr); else disp(A([6],:)); disp(lasterr); end 7(1) f=[]; for n=1:40

matlab程序设计第三章课后习题答案

1. p138 第6题在同一坐标轴中绘制下列两条曲线并标注两曲线交叉点。 >> t=0:0.01:pi; >> x1=t; >> y1=2*x1-0.5; >> x2=sin(3*t).*cos(t); >> y2=sin(3*t).*sin(t); >> plot(x1,y1,'r-',x2,y2,'g-') >> axis([-1,2,-1.5,1]) >> hold on >> s=solve('y=2*x-0.5','x=sin(3*t)*cos(t)','y=sin(3*t)*sin(t)'); >> plot(double(s.x),double(s.y),'*'); 截图:

p366 第4题绘制极坐标曲线,并分析对曲线形状的影响。 function [ output_args ] = Untitled2( input_args ) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here theta=0:0.01:2*pi; a=input('请输入a的值:'); b=input('请输入b的值:'); n=input('请输入n的值:'); rho=a*sin(b+n*theta); polar(theta,rho,'k'); end 下面以a=1,b=1,n=1的极坐标图形为基础来分析a、b、n的影响。

对a的值进行改变:对比发现a只影响半径值的整倍变化 对b的值进行改变:对比发现b的值使这个圆转换了一定的角度

对n的值进行改变:对比发现当n>=2时有如下规律 1、当n为整数时,图形变为2n个花瓣状的图形 2、当n为奇数时,图形变为n个花瓣状的图形 分别让n为2、3、4、5

MATLAB程序设计教程(第二版)第三章实验报告下载

大学社区网收集整理https://www.sodocs.net/doc/094730323.html, 评分 日期湖南商学院北津学院实验报告 课程名称MATLAB科学计算编程语言 实验名称MATLAB程序设计 专业班级信科1121班 姓名xxx 学号xxx 实验日期2012年11月5日 2012—2013学年度第一学期 一、实验目的 1.掌握利用if语句、switch语句实现选择结构的方法。 2.掌握利用for语句、while语句实现循环结构的方法。 3.熟悉利用向量运算来代替循环操作的方法并理解MATLAB程序设计的特点 4.掌握定义和调用MATLAB函数的方法。

二、实验环境 系统windows7旗舰版 处理器Intel(R)Core(TM)i7-3610M CPU @ 2.30GHz 安装内存 4.00GB (3.07GB 可用)系统类型64位操作系统运行环境 MATLAB 5.3 三、实验基本原理 利用上课所学知识解决以下问题: 1.从键盘输入一个3位数的整数,将它反向输出。如输入639,输出936。 2.输入一个百分制成绩,要求输出成绩等级A、B、C、D、E。其中90~100分为A,80~89分为B,70~79分为C,60~69分为D,60分以下为E。 要求: (1)分别用if 语句和switch 语句实现。 (2)输入百分制成绩后要判断该成绩的合理性,对不合理的成绩应输出出错信息。3.输入20个数,求其中最大数和最小数。要求分别用循环结构和调用MATLAB 的max 函数、min 函数来实现。 4.23.0ln )3.0sin(23.03.0a a e e y a a +++?=?,当a 取-3.0、-2.9、-2.8、…、2.8、2.9、 3.0时,求各点的函数值。要求分别用顺序结构和循环结构实现。 5.当n 分别取100、1000、10000时,求下列各式的值: (1)) 6...(n 1...31211122222π=+++++(2)) 2...()12)(12()2)(2(...756653443122π =??? ?????+?????????××????????××????????××n n n n 要求分别用循环结构和向量运算(使用sum 函数)来实现。 6.建立5×6矩阵,要求输出矩阵第n 行元素。当n 值超过矩阵的行数时,自动转为输出矩阵最后一行元素,并给出出错信息。 7已知,o999 ) 20()30()40(f f f y += (1)当)5ln(10)(2 ++=n n n f 时,y 的值是多小。 (2)当+×+×+×=433221)(n f …+)1(+×n n 时,y 的值是多小。 8.先用函数的递归调用定义一个函数文件求 ∑=n i m i 1,然后调用该函数文件求

实验3MATLAB程序设计

1,编写M 函数实现求一个数是否为素数,再编写一主程序(脚本文件),要求通过键盘输入一个整数,然后调用判断素数函数,从而确定它是否素数。 x=input('请输入一个整数x:'); if myprime(x) disp('您输入的整数x是一个素数。') else disp('您输入的数x不是一个素数。') end function y=myprime(x) y=1; for i=2:fix(sqrt(x)) if mod(x,i)==0 y=0; end end 2,编写M 函数统计一数值中零的个数,然后编写脚本文件,实现统计从1—2007 中零的总个数。 function num=number0(a) %统计十进制数值中0的个数 sa=num2str(a);%将数值装化为字符串 num=length(find(sa=='0'));% ));%求取字符串中'0’的个数 y=0;

for a=1:2006 num=number0(a); y=num+y; end disp(y) 504 3,编写程序计算x∈[-3,3],字长0.01:并画出曲线x = -3:0.01:3; y=zeros(size(x)); for i = 1:length(x) if -3<= x(i)& x(i)<=-1 y(i)=(-x(i).^2-4*x(i)-3)/ 2; elseif -1<= x(i) & x(i)<=1 y(i)=-x(i).^2+1; elseif 1<=x(:,i)<=3 y(i)=(-x(i).^2+4*x(i)-3)/2; end end plot(x,y) -3-2-10123

刘卫国版MATLAB程序设计与应用课后实验六八九

实验六 高层绘图操作 %第一题: 程序代码如下: x=linspace(0,2*pi,101); y=(0.5+3*sin(x)./(1+x.^2)).*cos(x); plot(x,y) 01234567 -1 -0.5 0.5 1 1.5 %第二题: %(1) 程序代码如下: x=linspace(-2*pi,2*pi,100); y1=x.^2; y2=cos(2*x); y3=y1.*y2; plot(x,y1,'b-',x,y2,'r:',x,y3,'y--'); text(4,16,'\leftarrow y1=x^2'); text(6*pi/4,-1,'\downarrow y2=cos(2*x)'); text(-1.5*pi,-2.25*pi*pi,'\uparrow y3=y1*y2');

-8 -6 -4 -2 2 4 6 8 -30-20 -10 10 20 30 40 %(2) 程序代码如下: x=linspace(-2*pi,2*pi,100); y1=x.^2; y2=cos(2*x); y3=y1.*y2; subplot(1,3,1);%分区 plot(x,y1); title('y1=x^2');%设置标题 subplot(1,3,2); plot(x,y2); title('y2=cos(2*x)'); subplot(1,3,3); plot(x,y3); title('y3=x^2*cos(2*x)');

-10 10 0510 15202530 35 40y1=x 2 -10 10 -1-0.8 -0.6 -0.4-0.200.20.4 0.6 0.8 1y2=cos(2*x) -10 10 -30-20 -10 10 20 30 40 y3=x 2*cos(2*x) %(3) 程序代码如下: x=linspace(-2*pi,2*pi,20); y1=x.^2; subplot(2,2,1);%分区 bar(x,y1); title('y1=x^2的条形图');%设置标题 subplot(2,2,2); stairs(x,y1); title('y1=x^2的阶梯图'); subplot(2,2,3); stem(x,y1); title('y1=x^2的杆图'); subplot(2,2,4); fill(x,y1,'r');%如果少了'r'则会出错 title('y1=x^2的填充图'); %其他的函数照样做。

实验二--MATLAB程序的设计(含实验报告)

实验二 MATLAB 程序设计 一、 实验目的 1.掌握利用if 语句实现选择结构的方法。 2.掌握利用switch 语句实现多分支选择结构的方法。 3.掌握利用for 语句实现循环结构的方法。 4.掌握利用while 语句实现循环结构的方法。 5.掌握MATLAB 函数的编写及调试方法。 二、 实验的设备及条件 计算机一台(带有MATLAB7.0以上的软件环境)。 M 文件的编写: 启动MATLAB 后,点击File|New|M-File ,启动MATLAB 的程序编辑及调试器(Editor/Debugger ),编辑以下程序,点击File|Save 保存程序,注意文件名最好用英文字符。点击Debug|Run 运行程序,在命令窗口查看运行结果,程序如有错误则改正 三、 实验容 1.编写求解方程02=++c bx ax 的根的函数(这个方程不一定为一元二次方程,因c b a 、、的不同取值而定),这里应根据c b a 、、的不同取值分别处理,有输入参数提示,当0~,0,0===c b a 时应提示“为恒不等式!”。并输入几组典型值加以检验。 (提示:提示输入使用input 函数) 2.输入一个百分制成绩,要求输出成绩等级A+、A 、B 、C 、D 、E 。其中100分为A+,90分~99分为A ,80分~89分为B ,70分~79分为C ,60分~69分为D ,60分以下为E 。 要求:(1)用switch 语句实现。 (2)输入百分制成绩后要判断该成绩的合理性,对不合理的成绩应输出出错信息。 (提示:注意单元矩阵的用法) 3.数论中一个有趣的题目:任意一个正整数,若为偶数,则用2除之,若为奇数,则与3相乘再加上1。重复此过程,最终得到的结果为1。如: 2→1 3→10→5→16→8→4→2→1 6→3→10→5→16→8→4→2→1

MATLAB程序设计实验报告

MATLAB 程序设计实验报告 一、实验目的 1. 通过实验熟悉MATLAB 仿真软件的使用方法; 2. 掌握用MATLAB 对连续信号时域分析、频域分析和s 域分析的方法,利用绘图命令绘制出典型信号的波形,了解这些信号的基本特征; 3. 掌握用MATLAB 对离散信号时域分析、频域分析和z 域分析的方法,利用绘图命令绘制出典型信号的波形,了解这些信号的基本特征; 4. 通过绘制信号运算结果的波形,了解这些信号运算对信号所起的作用。 二、实验设备 1. 计算机 : 2. MATLAB R2007a 仿真软件 三、实验原理 对系统的时域分析 信号的时域运算包括信号的相加、相乘,信号的时域变换包括信号的平移、反折、倒相及信号的尺度变换。 (1)信号的相加和相乘:已知信号)(1t f 和)(2t f ,信号相加和相乘记为 )()(1t f t f =)(2t f +;)()(1 t f t f =)(2t f *。 (2)信号的微分和积分:对于连续时间信号,其微分运算是用diff 函数来完成的,其语句格式为:diff(function,’variable’,n),其中function 表示需要进行求导运算的信号,或者被赋值的符号表达式;variable 为求导运算的独立变量;n 为求导的阶数,默认值为求一阶导数。连续信号的积分运算用int 函数来完成,语句格式为:diff(function,’variable’,a,b),其中function 表示需要进行被积信号,或者被赋值的符号表达式;variable 为求导运算的独立变量;a,b 为积分上、下限,a 和b 省略时为求不定积分。 (3)信号的平移、翻转和尺度变换 信号的平移包含信号的左移与右移,信号的翻转包含信号的倒相与折叠,平移和翻转信号不会改变信号)(t f 的面积和能量。信号的尺度变换是对信号)(t f 在时间轴上的变化,可使信号压缩或扩展。)(at f 将原波形压缩a 倍,)/(a t f 将原波形扩大a 倍。 ¥ 对系统频率特性的分析

Matlab编程与应用习题和一些参考答案

Matlab 上机实验一、二 3.求下列联立方程的解???????=+-+-=-+=++-=--+4 1025695842475412743w z y x w z x w z y x w z y x >> a=[3 4 -7 -12;5 -7 4 2;1 0 8 -5;-6 5 -2 10]; >> b=[4;4;9;4]; >> c=a\b 4.设???? ??????------=81272956313841A ,??????????-----=793183262345B ,求C1=A*B’;C2=A’*B;C3=A.*B,并求上述所有方阵的逆阵。 >> A=[1 4 8 13;-3 6 -5 -9;2 -7 -12 -8]; >> B=[5 4 3 -2;6 -2 3 -8;-1 3 -9 7]; >> C1=A*B' >> C2=A'*B >> C3=A.*B >> inv(C1) >> inv(C2) >> inv(C3) 5.设 ?? ????++=)1(sin 35.0cos 2x x x y ,把x=0~2π间分为101点,画出以x 为横坐标,y 为纵坐标的曲线。 >> x=linspace(0,2*pi,101); >> y=cos(x)*(0.5+(1+x.^2)\3*sin(x)); >> plot(x,y,'r') 6.产生8×6阶的正态分布随机数矩阵R1, 求其各列的平均值和均方差。并求该矩阵全体数的平均值和均方差。 (mean var ) a=randn(8,6) mean(a) var(a) k=mean(a) k1=mean(k) i=ones(8,6) i1=i*k1 i2=a-i1 i3=i2.*i2 g=mean(i3) g2=mean(g)

Matlab程序设计教程(第二版)刘卫国课后参考答案解析

第二章 1 求下列表达式的值。 (1) w=sqrt(2)*(1+0.34245*10^(-6)) (2) a=3.5; b=5; c=-9.8; x=(2*pi*a+(b+c)/(pi+a*b*c)-exp(2))/tan(b+c)+a (3) a=3.32; b=-7.9; y=2*pi*a^(2)*[(1-pi/4)*b-(0.8333-pi/4)*a] (4) t=[2,1-3*i;5,-0.65]; z=1/2*exp(2*t)*log(t+sqrt(1+t^(2))) 2 求下列表达式 A=[-1,5,-4;0,7,8;3,61,7]; B=[8,3,-1;2,5,3;-3,2,0]; (1) A+6*B A^2-B+eye (2) A*B A.*B B.*A (3) A/B B\A (4) [A,B] [A([1,3],:);B^2] 3 根据已知,完成下列操作 (1) A=[23,10,-0.778,0;41,-45,65,5;32,5,0,32;6,-9.54,54,3.14]; K=find(A>10&A<25); A(K) (2) A=[23,10,-0.778,0;41,-45,65,5;32,5,0,32;6,-9.54,54,3.14]; B=A(1:3,:) C=A(:,1:2) D=A(2:4,3:4) E=B*C (3) E

f=input('输入一个数:','s'); f(end :-1:1) 2 用if语句 score=input('请输入成绩:'); if score>=90&&score<=100 disp('A'); elseif score>=80&&score<=89 disp('B'); elseif score>=70&&score<=79 disp('C'); elseif score>=60&&score<=69; disp('D'); elseif score<60&&score>=0; disp('E'); else disp('出错'); end 用switch语句 score=input('请输入成绩:'); switch fix(score/10) case {9,10} disp('A'); case {8} disp('B'); case {7} disp('C'); case {6} disp('D'); case {0,1,2,3,4,5} disp('E'); otherwise disp('出错'); end 第四章1题 1) X=0:10; Y=x-x.^3/6; P lot(x,y)

Matlab程序设计实验报告

实验七Matlab 程序设计 实验目的: 1、掌握建立和执行M 文件的方法; 2、掌握实现选择结构的方法; 3、掌握实现循环结构的方法。 实验内容: 1. 编写用 5 次多项式拟合函数y=sin(x), x [0, 2 ]的脚本M 文件,要求绘图观察拟合的效果。 function shiyan1 x=0:0.5:2*pi y=sin(x) p=polyfit(x,y,5) x1=0:0.2:2*pi y1=polyval(p,x1) plot(x,y, 'b' ,x1,y1, '*r' x =

Columns 1 through 9 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 Columns 10 through 13 4.5000 5.0000 5.5000 6.0000 y = Columns 1 through 9 0 0.4794 0.8415 0.9975 0.9093 0.5985 0.1411 -0.3508 -0.7568 Columns 10 through 13 -0.9775 -0.9589 -0.7055 -0.2794 p = -0.0056 0.0881 -0.3967 0.2671 0.8902 0.0029 x1 = Columns 1 through 10 0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000 Columns 11 through 20

2. 2.2000 2.4000 2.6000 2.8000 3.0000 3.2000 3.4000 3.6000 1.8001 Columns 21 through 30 4.0 4.2000 4.4000 4.6000 4.8000 5.0000 5.2000 5.4000 5.6000 5.8000 Columns 31 through 32 6.0 6.2000 y1 = Columns 1 through 10 0.29 0.1886 0.3786 0.5585 0.7172 0.8461 0.9391 0.9926 1.0048 0.9761 Columns 11 through 20 0.9083 0.8048 0.6701 0.5098 0.3301 0.1381 -0.0590 -0.2538 -0.4389 -0.6073 Columns 21 through 30 -0.7524 -0.8685 -0.9505 -0.9949 -0.9991 -0.9626 -0.8863 -0.7732 -0.6288 -0.4606 Columns 31 through 32

matlab程序设计作业

Matlab程序设计作业 姓名: 学号: 专业:

? MATLAB 程序设计》作业 1、考虑如下x-y 一组实验数据: x=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10] y 二[1.2, 3, 4, 4, 5, 4.7, 5, 5.2, 6, 7.2] 分别绘出plot 的原始数据、一次拟合曲线和三次拟合曲线,给出 原始曲线 MATLAB 代码和运行结果。 7 6 5 4 3 2 2 3 4 5 6 7 8 9 10

7 6.5 6 5.5 5 4.5 4 3.5 3 2.5 10 一次拟合 三次拟合

x=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; y=[1.2, 3, 4, 4, 5, 4.7, 5, 5.2, 6, 7.2]; figure; plot(x,y) p1=polyfit(x,y,1); y1=polyval(p1,x); figure; plot(x,y1) p2=polyfit(x,y,3); y2=polyval(p2,x); figure; plot(x,y2) 2、在[0, 3n区间,绘制y二Sin(x)曲线(要求消去负半波,即(n 2n)区间内的函数值置零),求出曲线y 的平均值,以及y 的最大值及其最大值的位置。给出执行代码和运行结果。 x=0:pi/1000:3*pi; y=Sin(x); y1=(y>=0).*y; %消去负半波figure(1); plot(x,y1, 'b' ); a=mean(y1) %求出y1 的平均值 b=max(y1) %求出y1 的最大值b, 以及最大值在矩阵中的位置; d=x(find(y1==b)) >> ex1 a = 0.4243 b = 1 d = 1.5708 7.8540 >>

D实验五 M文件和MATLAB程序设计

实验五M文件和MATLAB程序设计 一、实验目的 matlab作为一种高级计算机语言,不仅可以命令行方式完成操作,也具有数据结构、控制流、输入输出等能力,本次实验通过熟悉和掌握m文件的建立与使用方法,以及函数与控制程序流程语句的使用,使学生具备一定的编程和程序调试能力。 1.掌握M文件的使用方法。 2.掌握if语句和switch语句的使用 3. 掌握循环语句的使用 4. 通过练习理解MATLAB编程方法。 二、实验原理 1.m文件 用matlab语言编写的程序,称为m文件。M文件根据调用方式的不同分为两类,命令文件(Script file)和函数文件(Function file)。区别? 2.程序控制结构 1)顺序结构 2)选择结构 (1)if语句a) 单分支if语句b) 双分支if语句c) 多分支if语句 (2)switch 语句 (3)try语句 3)循环结构 (1)for 语句 (2)while语句 (3)break语句、continue语句、return使用,区别? 3.函数文件 function 输出形参表=函数名(输入形参表) 注释说明部分 函数体语句 注意事项? 三、实验要求 1.首先上机练习PPT中各种流程控制语句的有关实例。 2.然后上机练习下面的实验习题。 四、实验习题

1.数论中一个有趣的题目:任意一个正整数,若为偶数,则用2除之,若为奇数,则与3相乘再加上1。重复此过程,最终得到的结果为1。如: 2→1 3→10→5→16→8→4→2→1 6→3→10→5→16→8→4→2→1 运行下面的程序,按程序提示输入n=1,2,3,5,7,8,9等数来验证这一结论。 %classic "3n+1" problem from number theory. while 1 n=input('Enter n,negative quits:'); if n<=0 break end a=n; while n>1 if rem(n,2)==0 n=n/2; else n=3*n+1; end a=[a,n]; end a end Enter n,negative quits: 2. 编程求满足∑=>m i i 11000020的最小m 值。 a=0; i=1; while (a<100000) a=a+pow2(i); i=i+1; end m=i-1 3. 编写一个函数,计算下面函数的值,给出x 的值,调用该函数后,返回y 的值。 function [y]=myfun1(x)

MATLAB程序设计教程课后答案.doc

MATLAB第二版课后答案unit3-8 unit3 实验指导 1、 n=input('请输入一个三位数:'); a=fix(n/100); b=fix((n-a*100)/10); c=n-a*100-b*10; d=c*100+b*10+a 2( 1) n=input(' 请输入成绩 '); switch n case num2cell(90:100) p='A'; case num2cell(80:89) p='B'; case num2cell(70:79) p='C'; case num2cell(60:69) p='D'; otherwise p='E'; end price=p (2) n=input(' 请输入成绩 '); if n>=90&n<=100 p='A'; elseif n>=80&n<=89 p='B'; elseif n>=70&n<=79 p='C'; elseif n>=60&n<=69 p='D'; else p='E'; end price=p (3) try n; catch price='erroe' end 3 n=[1,5,56,4,3,476,45,6,3,76,45,6,4,3,6,4,23,76,908,6];

a=n(1) ; b=n(1); for m=2:20 if n(m)>a a=n(m); elseif n(m)=0 disp(A([n],:)); elseif n<0 disp(lasterr); else disp(A([6],:)); disp(lasterr); end 7( 1) f=[];

实验二 MATLAB程序设计

实验二 MATLAB 程序设计 一、 实验目的 1.掌握利用if 语句实现选择结构的方法。 2.掌握利用switch 语句实现多分支选择结构的方法。 3.掌握利用for 语句实现循环结构的方法。 4.掌握利用while 语句实现循环结构的方法。 5.掌握MATLAB 函数的编写及调试方法。 二、 实验的设备及条件 计算机一台(带有以上的软件环境)。 M 文件的编写: 启动MATLAB 后,点击File|New|M-File ,启动MATLAB 的程序编辑及调试器(Editor/Debugger ),编辑以下程序,点击File|Save 保存程序,注意文件名最好用英文字符。点击Debug|Run 运行程序,在命令窗口查看运行结果,程序如有错误则改正 三、 实验内容 1.编写求解方程02=++c bx ax 的根的函数(这个方程不一定为一元二次方程,因c b a 、、的不同取值而定),这里应根据c b a 、、的不同取值分别处理,有输入参数提示,当0~,0,0===c b a 时应提示“为恒不等式!”。并输入几组典型值加以检验。 (提示:提示输入使用input 函数) 2.输入一个百分制成绩,要求输出成绩等级A+、A 、B 、C 、D 、E 。其中100分为A+,90分~99分为A ,80分~89分为B ,70分~79分为C ,60分~69分为D ,60分以下为E 。 要求:(1)用switch 语句实现。 (2)输入百分制成绩后要判断该成绩的合理性,对不合理的成绩应输出出错信息。 (提示:注意单元矩阵的用法) 3.数论中一个有趣的题目:任意一个正整数,若为偶数,则用2除之,若为奇数,则与3相乘再加上1。重复此过程,最终得到的结果为1。如: 21 21 421 运行下面的程序,按程序提示输入n=1,2,3,5,7等数来验证这一结论。 请为关键的Matlab 语句填写上相关注释,说明其含义或功能。 4. y 的值。 5. (Root Mean Square)的计算公式为:

MATLAB程序设计与应用(第二版)实验参考答案

MATLAB程序设计与应用(第二版)实验参考答案 %实验一MATLAB运算基础 %第一题ftp://192.168.0.143/ %(1) z1=2*sin(85*pi/180)/(1+exp(2))ftp://192.168.0.143/ %(2) x=[2,1+2i;-0.45,5]; z2=0.5*log(x+sqrt(1+x.^2))ftp://192.168.0.143/ %(3) a=-3.0:0.1:3.0;ftp://192.168.0.143/ z3=(exp(0.3*a)-exp(-0.3*a))/2.*sin(a+0.3)+log((0.3+a)/2) %(4) t=0:0.5:2.5; z4=t.^2.*(t>=0&t<1)+(t.^2-1).*(t>=1&t<2)+(t.^2-2*t+1).*(t>=2&t<3) %第二题 A=[12 34 -4;34 7 87;3 65 7]; B=[1 3 -1;2 0 3;3 -2 7]; A+6*B A-B+eye(size(A)) A*B A.*B A^3 A.^3 A/B B\A [A,B] [A([1,3],:);B^2] %第三题 A=[1 2 3 4 5;6 7 8 9 10;11 12 13 14 15;16 17 18 19 20;21 22 23 24 25] B=[3 0 16;17 -6 9;0 23 -4;9 7 0;4 13 11] C=A*B F=size(C) D=C(F(1)-2:F(1),F(2)-1:F(2)) whos %第四题 %(1): A=100:999; B=rem(A,21); C=length(find(B==0)) %(2): A='lsdhKSDLKklsdkl';

Matlab程序设计(2016大作业)

Matlab程序设计 课程大作业 题目名称:_________________________________ 班级:_________________________________ 姓名:_________________________________ 学号:_________________________________ 课程教师:温海骏 学期:2015-2016学年第2学期 完成时间: MATLAB优化应用 §1 线性规划模型 一、线性规划问题: 问题1:生产计划问题 假设某厂计划生产甲、乙两种产品,现库存主要材料有A类3600公斤,B类2000公斤,C类3000公斤。每件甲产品需用材料A类9公斤,B类4公斤,C类3公斤。每件乙产品,需用材料A类4公斤,B类5公斤,C类10公斤。甲单位产品的利润70元,乙单位产品的利润120元。问如何安排生产,才能使该厂所获的利润最大。 问题2:投资问题 某公司有一批资金用于4个工程项目的投资,其投资各项目时所得的净收益(投入资金百分比)如下表:工程项目收益表 工程项目 A B C D 收益(%) 15 10

12 由于某种原因,决定用于项目A的投资不大于其他各项投资之和而用于项目B和C的投资要大于项目D的投资。试确定该公司收益最大的投资分配方案。 问题3:运输问题 有A、B、C三个食品加工厂,负责供给甲、乙、丙、丁四个市场。三个厂每天生产食品箱数上限如下表: 工厂 A B C 生产数 60 40 50 四个市场每天的需求量如下表: 市场 甲 乙 丙 丁 需求量 20 35 33 34 从各厂运到各市场的运输费(元/每箱)由下表给出: 收点 发点 市场 甲 乙 丙 丁 工 厂 A 2 1 3 2 B