Matlab源程序
matlab实现数值计算功能源程序(个人整理)

matlab数值计算功能1,基础运算(1)多项式的创建与表达将多项式(x-6)(x-3)(x-8)表示为系数形式a=[6 3 8] % 写成根矢量pa=poly(a)% 求出系数矢量ppa=poly2sym(pa,'x') % 表示成符号形式ezplot(ppa,[-50,50])求3介方阵A的特征多项式a=[6 2 4;7 5 6;1 3 6 ];pa=poly(a)% 写出系数矢量ppa=poly2sym(pa) %表示成符号形式ezplot(ppa,[-50,50]) % 绘图求x^3-6x^2-72x-27的根。
a=[1,-6,-72,-85]; % 写出多项式系数矢量r=roots(a) % 求多项式的根(2)多项式的乘除运算c=conv(a,b) %乘法[q,r]=deconv(c,a)% 除法求a(s)=s^2+2s+3乘以b(s)=4s^2+5s+6的乘积a=[1 2 3]b=[4 5 6] % 写出系数矢量c=conv(a,b)c=poly2sym(c,'s') % 写成符号形式的多项式展开(s^2+2s+2)(s+4)(s+1)并验证结果被(s+4),(s+3)除后的结果。
c=conv([1,2,2],conv([1,4],[1,1]))cs=poly2sym(c,'s')c=[1 7 16 18 8][q1,r1]=deconv(c,[1,4])[q2,r2]=deconv(c,[1,3])cc=conv(q2,[1,3])test=((c-r2)==cc)其他常用的多项式运算命令pa=polyval(p,s) % 按数组规则计算给定s时多项式的值pm=polyvalm(p,s)% 按矩阵规则计算给定s时多项式的值[r,p,k]=residue(b,a) % 部分分式展开,b,a分别是分子,分母多项式系数矢量。
r,p,k分别是留数,极点和值项矢量。
OFDM_matlab源程序总结

比较完整的OFDM仿真,供大家学习下载。
、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、%%一个相对完整的OFDM通信系统的仿真设计,包括编码,调制,IFFT,%上下变频,高斯信道建模,FFT,PAPR抑制,各种同步,解调和解码等模%块,并统括系统性能的仿真验证了系统设计的可靠性。
clear allclose allclc%++++++++++++++++++++++++++全局变量++++++++++++++++++++++++++++++% seq_num 表示当前帧是第几帧% count_dds_up 上变频处的控制字的累加% count_dds_down 下变频处的控制字的累加(整整)% count_dds_down_tmp 下变频处的控制字的累加(小数)% dingshi 定时同步的定位% m_syn 记录定时同步中的自相关平台global seq_numglobal count_dds_upglobal count_dds_downglobal count_dds_down_tmpglobal dingshiglobal m_syn%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++% SNR_Pre 设定用于仿真的信噪比的初值% interval_SNR 设定用于仿真的信噪比间隔% frame_num 每一个信噪比下仿真的数据帧数% err_int_final 用于计算每帧出现的误比特数% fwc_down 设定的接收机初始载波频率控制字% fre_offset 设定接收机初始载波频率偏移调整量(单位为Hz)% k0 每次进入卷积编码器的信息比特数% G 卷积编码的生成矩阵SNR_Pre=-5;interval_SNR=1;for SNR_System=SNR_Pre:interval_SNR:5frame_num=152;dingshi=250;err_int_final=0;fwc_down=;fre_offset=0;k0=1;G=[1 0 1 1 0 1 1;1 1 1 1 0 0 1 ];disp('--------------start-------------------');for seq_num=1:frame_num, %frame_num 帧数%+++++++++++++++++++++++以下为输入数据部分+++++++++++++++++++++++datain=randint(1,90);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++++以下为信道卷积编码部分+++++++++++++++++++++ encodeDATA=cnv_encd(G,k0,datain);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++信道交织编码+++++++++++++++++++++++++interlacedata=interlacecode(encodeDATA,8,24);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++以下为QPSK调制部分+++++++++++++++++++++QPSKdata=qpsk(interlacedata);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++生成训练序列+++++++++++++++++++++++++if seq_num<3trainsp_temp=seq_train();end%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++++++++++插入导频++++++++++++++++++++++++++++PILOT=(1+j);m_QPSKdata=QPSKdata;data2fft_temp=[m_QPSKdata(1:8),PILOT,m_QPSKdata(9:16),PILOT,m_QPSKdata(17:24),P ILOT,m_QPSKdata(25:32),PILOT,m_QPSKdata(33:40),PILOT,m_QPSKdata(41:48),m_QPSKda ta(49:56),PILOT,m_QPSKdata(57:64),PILOT,m_QPSKdata(65:72),PILOT,m_QPSKdata(73:8 0),PILOT,m_QPSKdata(81:88),PILOT,m_QPSKdata(89:end)];%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++trainsp_temp2=[trainsp_temp,zeros(1,128)];trainsp=[trainsp_temp2(65:256),trainsp_temp2(1:64)];%++++++++++++++++++++++++++降PAPR矩阵变换++++++++++++++++++++++++matix_data=nyquistimp_PS();matrix_mult=data2fft_temp*matix_data;%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++data2fft2=[matrix_mult(65:128),zeros(1,128),matrix_mult(1:64)];%++++++++++++++++++++++++++++ifft运算+++++++++++++++++++++++++++if seq_num==1ifftin=trainsp;elseif seq_num==2ifftin=trainsp;elseifftin=data2fft2;endIFFTdata=fft_my(conj(ifftin)/256);IFFTdata=conj(IFFTdata);% figure% plot(real(IFFTdata))% xlabel('realIFFTdata')% figure% plot(imag(IFFTdata))% xlabel('imagIFFTdata')%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++以下为插入循环前后缀,2倍升采样+++++++++++++++data2fir=add_CYC_upsample(IFFTdata,2);% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++% +++++++++++++++++++++++++fir低通滤波++++++++++++++++++++++++++guiyi_a=[];%抽样截止频率为128kHZ,通带截止频率为20kHZ,阻带截止频率为40kHZ,带内纹波动小于1dB,带外衰减100dBtxFIRdatai=filter(guiyi_a,1,real(data2fir));txFIRdataq=filter(guiyi_a,1,imag(data2fir));% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++发射机cic滤波++++++++++++++++++++++++++CICidatai=cic_inter(txFIRdatai,20);CICidataq=cic_inter(txFIRdataq,20);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++++上变频++++++++++++++++++++++++++++fwc_up=16; %控制字可以选择DUCdata=up_convert_ofdm(fwc_up,CICidatai,CICidataq);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++++++++高斯白噪声信道++++++++++++++++++++++++[DUCdata,datamax]=guiyi_DUCdata(DUCdata);awgn_data=awgn(DUCdata,SNR_System);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%*****************************接受机*****************************%+++++++++++++++++++++++++++++下变频+++++++++++++++++++++++++++++DUCdata_tmp=awgn_data;fwc_down=fwc_down+(fre_offset*128/2560000);r_fre_offset=2560000*((fwc_down-fwc_up)/128);[DDCdatai,DDCdataq]=down_convert_ofdm(fwc_down,DUCdata_tmp);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++++++++接收机cic滤波+++++++++++++++++++++++++CICddatai=cic_deci(DDCdatai,40,40);CICddataq=cic_deci(DDCdataq,40,40);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++++++++++fir低通滤波++++++++++++++++++++++++guiyi_b=[ ];%抽样截止频率为64kHZ,通带截止频率为20kHZ,阻带截止频率为30kHZ,带内纹波动小于1dB,带外衰减60dBrxFIRdatai=filter(guiyi_b,1,CICddatai);rxFIRdataq=filter(guiyi_b,1,CICddataq);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++++++量化++++++++++++++++++++++++++++q_rxFIRdatai=sign(rxFIRdatai);q_rxFIRdataq=sign(rxFIRdataq);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++++++++++定时同步检测++++++++++++++++++++++++if seq_num<3time_syn(q_rxFIRdatai,q_rxFIRdataq);end%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++++频率同步+++++++++++++++++++++++++++fre_offset=fre_syn(rxFIRdatai,rxFIRdataq);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++++fft运算+++++++++++++++++++++++++++if seq_num>2seq_num-2fftw=32+dingshi;rxFIRdata_syn=rxFIRdatai(fftw:fftw+255)+j*rxFIRdataq(fftw:fftw+255); FFTdata=fft_my(rxFIRdata_syn);%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++++++++降PAPR逆矩阵变换++++++++++++++++++++++ fftdata_reg=[FFTdata(193:256),FFTdata(1:64)];dematrix_data=fftdata_reg*pinv(matix_data);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++++++++++相位补偿+++++++++++++++++++++++++++rx_qpsk_din_th=phase_comp(dematrix_data);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++QPSK解调部分++++++++++++++++++++++++% figure% plot(rx_qpsk_din_th,'.')% xlabel('星座图')datatemp4=deqpsk(rx_qpsk_din_th);datatemp4=sign(datatemp4);for m=1:192if datatemp4(m)==-1datatemp4(m)=1;elseif datatemp4(m)==1datatemp4(m)=0;endend%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++信道解交织++++++++++++++++++++++++++++ interdout=interlacedecode(datatemp4,8,24);%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%++++++++++++++++++++以下为viterbi译码部分++++++++++++++++++++++ decodeDATA=viterbi(G,k0,interdout);%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++误比特统计++++++++++++++++++++++++++++err_final=sum(abs(decodeDATA-datain))err_int_final=err_int_final+err_finalendenddisp('------------------------------------------------------');SNR_Systemerr_rate_final((SNR_System-SNR_Pre)./interval_SNR+1)=err_int_final/(90*(frame_n um-2))disp('------------------------------------------------------');enddisp('--------------end-------------------');SNR_System=SNR_Pre:interval_SNR:5;figuresemilogy(SNR_System,err_rate_final,'b-*');xlabel('信噪比/dB')ylabel('误码率')axis([-5,5,0,1])grid on%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%************************beginning of file*****************************%%卷积码编码程序function output=cnv_encd(G,k0,input)% cnv_encd(G,k0,input),k0 是每一时钟周期输入编码器的 bit 数,% G 是决定输入序列的生成矩阵,它有 n0 行 L*k0 列 n0 是输出 bit 数,% 参数 n0 和 L 由生成矩阵 G 导出,L 是约束长度。
数学建模2020c题matlab源程序

数学建模2020c题matlab源程序摘要:1.介绍数学建模2020c 题的背景和意义2.概述Matlab 在解决该题中的应用3.分析题目要求的解决方案4.提供Matlab 源程序的实现方法5.总结Matlab 在解决该题中的优势和意义正文:一、介绍数学建模2020c 题的背景和意义数学建模2020c 题是一道具有挑战性的数学问题,它涉及到许多数学领域的知识,如微积分、线性代数、概率论等。
这道题目旨在考查学生运用数学知识解决实际问题的能力,通过构建数学模型,找到问题的解决方案。
在解决这道题目的过程中,Matlab 作为一种强大的数学软件,可以提供很大的帮助。
二、概述Matlab 在解决该题中的应用Matlab 是一种功能强大的数学软件,它可以用于解决各种数学问题,包括数学建模问题。
在解决数学建模2020c 题时,Matlab 可以帮助学生进行数据分析、绘制图像、求解方程等。
通过使用Matlab,学生可以更轻松地构建数学模型,找到问题的解决方案。
三、分析题目要求的解决方案为了解决数学建模2020c 题,首先需要对题目进行仔细的分析。
题目要求我们找到一个数学模型,可以描述问题的解决方案。
通过分析题目,我们可以发现,这道题目需要我们运用微积分、线性代数、概率论等知识,构建一个数学模型,求解方程,找到问题的解决方案。
四、提供Matlab 源程序的实现方法在解决数学建模2020c 题时,Matlab 可以提供很大的帮助。
通过使用Matlab,我们可以轻松地构建数学模型,求解方程,找到问题的解决方案。
下面是一个简单的Matlab 源程序,可以实现题目要求的解决方案。
```Matlab% 初始化参数a = 1;b = 2;c = 3;x0 = 0;y0 = 0;t0 = 0;% 求解微分方程df = @(x, t) [a * x + b * y; c * x - a * y];df_init = [x0; y0];[x, y] = ode45(@df, [0, t0], df_init);% 绘制图像figure;plot3(x, y, t);xlabel("x");ylabel("y");zlabel("t");title("Solution of the problem");```五、总结Matlab 在解决该题中的优势和意义通过使用Matlab,我们可以更轻松地解决数学建模2020c 题。
调用Matlab对矩阵做奇异值分解的程序源代码

调用Matlab对矩阵做奇异值分解的程序源代码为计算针对多目标节点集的相对连通系数需要对矩阵进行奇异值分解,为此可直接调用Matlab里的对应矩阵计算模块。
从VC++调用Matlab进行计算有两种方法:一种是通过Matlab引擎直接调用Matlab,另一种是将Matlab函数编译为dll在VC中使用。
这里采用第二种方法,通过Matlab引擎调用矩阵计算功能。
其实现代码如下。
void Connectivitys::svd_Matlab(void * destMatrixP){int i,j;double * vP=NULL,* sP=NULL;mxArray *Xin,*Vout,*Sout;Engine * ep;//显式加载Matlab动态链接库,定义其中的函数HINSTANCE hInst_libeng=LoadLibrary(“libeng.dll”);HINSTANCE hInst_libmat=LoadLibrary(“libmat.dll”);HINSTANCE hInst_libmx=LoadLibrary(“libmx.dll”);typedef Engine *(*_EngOpen)(const char *startcmd);_EngOpen engOpen=(_EngOpen)GetProcAddress(hInst_libeng,“engOpen”);typedef int(*_EngClose)(Engine*ep);_EngClose engClose=(_EngClose)GetProcAddress(hInst_libeng,“engClose”);typedef mxArray *(*_EngGetArray)(Engine *ep,const char *name);_EngGetArray engGetArray=(_EngGetArray)GetProcAddress(hInst_libeng,“engGetArray”);typedef int(*_EngPutArray)(Engine *ep,const mxArray *ap);_EngPutArray engPutArray=(_EngPutArray)GetProcAddress(hInst_libeng,“engPutArray”);typedef int(*_EngEvalString)(Engine *ep,const char *string);_EngEvalString engEvalString=(_EngEvalString)GetProcAddress(hInst_libeng,“engEvalString”);typedef mxArray *(*_MxCreateDoubleMatrix)(int m,int n,mxComplexity flag);_MxCreateDoubleMatrix mxCreateDoubleMatrix=(_MxCreateDoubleMatrix)GetProcAddress(hInst_libmx,“mxCreateDoubleMatrix”);typedef void(*_MxDestroyArray)(mxArray *pa);_MxDestroyArray mxDestroyArray=(_MxDestroyArray)GetProcAddress(hInst_libmx,“mxDestroyArray”);typedef void(*_MxSetName)(mxArray*pa,const char *s);_MxSetName mxSetName=(_MxSetName)GetProcAddress(hInst_libmx,“mxSetName”);typedef double *(*_MxGetPr)(const mxArray *pa);_MxGetPr mxGetPr=(_MxGetPr)GetProcAddress(hInst_libmx,“mxGetPr”);//调用MATLAB引擎if(!(ep=engOpen(NULL))){AfxMessageBox(_T(“Can’t start MATLAB engine”));exit(-1);}//将平均偏差矩阵载入MatlabXin=mxCreateDoubleMatrix(destNodesNum,totalNodeNum,mxREAL);mxSetName(Xin,“x”);memcpy((char *)mxGetPr(Xin),(char *)destMatrixP,destNodesNum*totalNodeNum*sizeof(double));free(destMatrixP);engPutArray(ep,Xin);//对平均偏差矩阵进行SVD分解engEvalString(ep,“[u,s,v]=svd(x);”);//获取SVD分解结果Sout=engGetArray(ep,“s”);sP=mxGetPr(Sout);Vout=engGetArray(ep,“v”);vP=mxGetPr(Vout);//利用SVD分解结果进行多尺度分析//求多目标节点集相对连通系数RCCfloat totalSum=0.0,vi=0.0;for(j=0;j<destNodesNum;j++)totalSum+=*(sP+(destNodesNum+1)*j);for(i=0;i<totalNodeNum;i++){MultiNodesRCC[i]=0.0;for(j=0;j<destNodesNum;j++){vi=*(vP+totalNodeNum*j+i);if(vi<0)vi=0-vi;MultiNodesRCC[i]+=(*(sP+(destNodesNum+1)*j)/totalSum)*vi;}}//关闭MATLAB引擎,释放分配的资源engClose(ep);mxDestroyArray(Xin);mxDestroyArray(Sout);mxDestroyArray(Vout);FreeLibrary(hInst_libeng);FreeLibrary(hInst_libmat);FreeLibrary(hInst_libmx);}-全文完-。
数字下变频matlab源程序

数字下变频matlab源程序当涉及到数字信号处理中的频率变换,MATLAB是一个非常强大的工具。
在MATLAB中,可以使用不同的函数和工具箱来实现数字信号的频率变换。
下面是一个简单的示例,展示了如何在MATLAB中实现数字信号的频率变换。
matlab.% 生成输入信号。
fs = 1000; % 采样频率。
t = 0:1/fs:1-1/fs; % 时间向量。
f1 = 20; % 输入信号频率。
x = sin(2pif1t); % 输入信号。
% 进行频率变换。
f2 = 50; % 目标频率。
y = x.exp(1i2pif2t); % 频率变换。
% 绘制结果。
subplot(2,1,1);plot(t,real(x));title('原始信号');xlabel('时间');ylabel('幅度');subplot(2,1,2);plot(t,real(y));title('频率变换后的信号');xlabel('时间');ylabel('幅度');在这个示例中,我们首先生成了一个输入信号x,然后使用频率变换公式y = x.exp(1i2pif2t)将输入信号的频率变换到f2。
最后,我们绘制了原始信号和频率变换后的信号的波形图。
除了这个简单的示例之外,MATLAB还提供了许多内置的函数和工具箱,如fft, ifft, chirp, spectrogram等,可以用于数字信号的频率变换。
通过这些函数和工具箱,可以实现更复杂和高级的频率变换操作,比如滤波、混频、调制解调等。
总之,MATLAB是一个非常强大的工具,可以帮助你实现数字信号的频率变换。
希望这个简单的示例可以帮助你更好地理解在MATLAB中实现数字信号频率变换的方法。
因子分析MATLAB程序源代码

因子分析MATLAB程序源代码因子分析是一种统计分析方法,用于确定一组观测变量的潜在因子结构。
在MATLAB中,可以使用`factoran`函数进行因子分析。
下面是一个示例的MATLAB代码,用于执行因子分析,并使用经典的因子旋转方法-变简单结构法(varimax)进行因子旋转:```matlab%数据准备data = [1 2 3 4 5; 2 3 4 5 6; 3 4 5 6 7; 4 5 6 7 8; 5 6 7 8 9];%因子分析num_factors = 2; % 欲提取的因子个数method = 'varimax'; % 因子旋转方法[loadings, spec_var, T, stats] = factoran(data, num_factors, 'rotate', method);%输出结果disp('因子载荷矩阵:');disp(loadings);disp('特殊因子方差:');disp(spec_var);disp('变换矩阵:');disp(T);disp('统计信息:');disp(stats);```在这个示例中,我们首先准备了一个包含5个观测变量的数据矩阵`data`,然后指定欲提取的因子个数`num_factors`为2,并选择`varimax`作为因子旋转方法。
调用`factoran`函数进行因子分析后,会返回以下结果:- `loadings`:因子载荷矩阵,描述所有观测变量与因子之间的关系;- `spec_var`:特殊因子方差,每个观测变量独自解释的方差;-`T`:变换矩阵,将数据旋转到简单结构后的载荷矩阵;- `stats`:统计信息,包括共同度、特殊方差和解释总方差等。
使用`disp`函数将结果打印到控制台上。
以上就是一个简单的因子分析的MATLAB代码示例,你可以根据自己的需求修改和扩展代码。
数学建模2020c题matlab源程序

数学建模2020c题是一项涉及到数学、计算机等多个学科知识的综合性竞赛题目。
在这次比赛中,参赛者需要应用数学建模方法,解决实际问题,并编写相应的程序进行模拟和计算。
其中,matlab作为一种强大的数学软件工具,被广泛应用于数学建模的程序编写中。
一、问题分析在数学建模2020c题中,通常会给出一个实际问题,同时要求参赛者利用数学建模的方法进行分析和求解。
在编写matlab源程序时,需要遵循以下几个步骤:1. 问题理解与建模: 首先要对问题进行深入理解,明确问题的要求和限制条件,然后进行数学建模,将实际问题转化为数学模型。
2. 程序框架设计: 根据数学模型的建立,设计matlab程序的框架结构,确定程序的主要功能和模块划分。
3. 算法设计与实现: 针对数学模型所需的运算和计算过程,设计相应的算法,并在matlab中进行实现。
4. 程序测试与优化: 编写完程序后,需要进行充分的测试,验证程序的正确性和稳定性,并对程序进行优化,提高程序的运行效率。
二、matlab源程序编写在数学建模竞赛中,matlab源程序往往会涉及到各种数学方法和算法,如数值计算、优化算法、模拟仿真等。
下面以一个示例问题为例,介绍matlab源程序的编写过程。
示例题目:某市的交通拥堵问题问题描述:某市的交通拥堵问题日益严重,为了解决交通拥堵问题并提高交通效率,需要对该市的交通流进行合理规划和控制。
假设该市的道路网络和车辆流量情况已知,请设计一个算法,模拟和优化该市的交通流,并给出相应的交通规划建议。
1. 问题理解与建模根据问题描述,需要对该市的道路网络和车辆流量进行建模。
首先将道路网络抽象为一个图结构,节点表示道路交叉口,边表示道路连通关系;然后对车辆的流动进行建模,考虑车辆的速度、车流密度等因素。
2. 程序框架设计根据问题的数学模型,设计matlab程序的框架结构。
程序主要包括道路网络的表示与初始化、车辆流动模拟与优化算法、交通规划建议输出等模块。
magnify matlab源程序

magnify matlab源程序以下是一个简单的 MATLAB 源程序,用于实现图片的放大功能: % 以下代码实现了图片的放大功能% 输入的变量 im 是一个三维 RGB 图像矩阵,factor 是放大因子 % 输出的变量 im_scaled 是一个放大后的三维 RGB 图像矩阵function im_scaled = magnify(im, factor)% 读取输入的图片大小[height, width, ~] = size(im);% 通过放大因子计算新的图片大小new_height = round(height * factor);new_width = round(width * factor);% 通过双线性插值法计算新的像素值,并把结果存储到输出矩阵中 for i = 1:new_heightfor j = 1:new_widthy = y_ratio * (i-1) + 1;x = x_ratio * (j-1) + 1;y1 = floor(y);y2 = ceil(y);x1 = floor(x);x2 = ceil(x);if x1 == x2x2 = x1 + 1;endif y1 == y2y2 = y1 + 1;endf11 = double(im(y1, x1, :));f12 = double(im(y2, x1, :));f21 = double(im(y1, x2, :));f22 = double(im(y2, x2, :));dx = x - x1;dy = y - y1;im_scaled(i,j,:) = uint8((1-dx)*(1-dy)*f11 + dx*(1-dy)*f21 +(1-dx)*dy*f12 + dx*dy*f22);endend% 显示原始图片和放大后的图片subplot(1, 2, 1);imshow(im);title('原始图片');subplot(1, 2, 2);imshow(im_scaled);title('放大后的图片');end该程序的功能是将输入的 RGB 图片矩阵进行放大,输出结果是放大后的图片矩阵。
matlab ode78源程序

MATLAB是一种广泛用于科学计算和工程设计的高级编程语言和交互式环境。
其ODE78函数是MATLAB中用于求解常微分方程组的函数,能够高效地求解一般形式的刚性和非刚性常微分方程组。
本文将介绍ODE78函数的源程序以及其相关知识。
一、ODE78函数介绍1. ODE78函数是MATLAB中的一个常微分方程(ODE)求解器,可用于解决形如dy/dx=f(x,y)的一阶和高阶ODE。
2. ODE78函数利用Adams-Bashforth-Moulton方法和变步长算法实现了高效的ODE求解。
3. ODE78函数支持对ODE方程的自定义,并可通过设置选项来调整算法参数以获得更精确的求解结果。
二、ODE78函数源程序1. ODE78函数的源程序如下所示:function [T,Y] = ode78(odefun,tspan,y0,options)t0 = tspan(1); tf = tspan(length(tspan));hmax=0.1; tol=1.e-3;if nargin < 4 options=[]; endif ~isempty(options)if isfield(options,'hmax') hmax=options.hmax; endif isfield(options,'tol') tol=options.tol; endendh=0.1; y=y0(:); t = t0; T = t; Y = y'; iflag=1;m本人n loopwhile iflag > 0h0=h;if t+h > tfh = tf - t;endk1=feval(odefun,t,y);k2=feval(odefun,t+h/7,y+(h/7)*k1);k3=feval(odefun,t+h/6,y+(h/6)*k1+(h/6)*k2);k4=feval(odefun,t+(5*h/12),y+(5*h/12)*k1-(25*h/16)*k2+(25*h/16)*k3);k5=feval(odefun,t+(3*h/4),y+(3*h/4)*k1 +(3*h/4)*k2-(3*h/4)*k3+(9*h/4)*k4); k6=feval(odefun,t+h,y+(h)*k1-(3*h)*k2+(12*h)*k3-(12*h)*k4+(8*h)*k5);k2p=feval(odefun,t+(7*h/8),y+(7*h/8)*k1 + (7*h/8)*k2 -(21*h/8)*k3 +(35*h/8)*k4 -(35*h/8)*k5 +(7*h/8)*k6);k5p=feval(odefun,t+(3*h/4),y+(3*h/4)*k1 +(3*h/4)*k2-(3*h/4)*k3+(9*h/4)*k4 +(9*h/4)*k5-(3*h/4)*k6);ynew=y+h*(k1/11+28*k3/54+31*k4/54-2*k5/27); k5=k5-k5p; y5=tol/max(abs(y),ynew.*tol)+max(abs(y),ynew.*tol);if max(y5)-1<= .OKy=ynew; t=t+h; T=[T t]; Y=[Y; y']; endend2. 上文是ODE78函数的源程序,其中包括了主要的计算过程和参数设置。
matlab中bode函数源程序 -回复

matlab中bode函数源程序-回复如何使用MATLAB中的bode函数绘制频率响应曲线MATLAB是一种功能强大的数值计算与数据可视化软件,广泛应用于科学与工程领域。
在信号处理和控制系统设计中,频率响应曲线是一项重要的分析工具。
MATLAB中的bode函数可以方便地绘制频率响应曲线,并帮助工程师更好地了解系统的频域特性。
本文将介绍如何使用MATLAB中的bode函数来绘制频率响应曲线,让我们一步一步来进行解答。
第一步:了解bode函数的基本用法在开始使用bode函数之前,我们需要了解其基本用法。
bode函数是MATLAB中的一个内置函数,用于绘制线性时不变系统的频率响应曲线。
该函数的基本语法如下:[bode_mag, bode_phase, wout] = bode(sys)其中,sys是一个已定义的线性时不变系统,bode_mag是频率响应曲线的幅值,bode_phase是频率响应曲线的相位,wout是频率向量。
第二步:创建线性时不变系统在使用bode函数之前,我们需要创建一个线性时不变系统。
有多种方式可以创建一个线性时不变系统,例如传递函数(tf),状态空间模型(ss),零极点(zpk)等。
以下是一个通过传递函数创建线性时不变系统的例子:sys = tf(num,den)其中,num是传递函数的分子系数,den是传递函数的分母系数。
第三步:确定频率范围在使用bode函数绘制频率响应曲线之前,我们需要确定画图的频率范围。
可以通过设定一个频率向量来实现,例如:w = logspace(-3,3,1000)其中,logspace函数用于生成一个以10为底数的对数均匀分布的向量。
上述例子中,我们生成了一个从10^-3到10^3的频率向量,总共包含1000个点。
第四步:调用bode函数绘制频率响应曲线一旦我们创建了线性时不变系统和确定了频率范围,就可以调用bode 函数来绘制频率响应曲线了。
以下是一个完整的绘制频率响应曲线的例子:sys = tf([1],[1 2 1])w = logspace(-3,3,1000)[bode_mag, bode_phase, wout] = bode(sys,w)在这个例子中,我们创建了一个二阶系统1/(s^2 + 2s + 1),设定了频率范围为10^-3到10^3,并调用bode函数绘制频率响应曲线。
IIR滤波器matlab源程序

IIR滤波器matlab源程序(1)IIR一阶低通滤波器clear;fi=1;fs=10;Gc2=0.9;wc=2*pi*fi/fs;omegac=tan(wc/2);alpha=(sqrt(Gc2)/sqrt(1-Gc2))*omegac;a=(1-alpha)/(1+alpha);b=(1-a)/2;w=0:pi/300:pi;Hw2=alpha^2./(alpha^2+(tan(w/2)).^2);plot(w/pi,Hw2);grid;hold on;(2)一阶高通滤波器clear;fi=1;fs=10;Gc2=0.5;wc=2*pi*fi/fs;omegac=tan(wc/2);alpha=(sqrt(1-Gc2)/(sqrt(Gc2)))*omegac;a=(1-alpha)/(1+alpha);b=(1+a)/2;w=0:pi/300:pi;Hw2=(tan(w/2).^2)./(alpha^2+(tan(w/2)).^2); plot(w/pi,Hw2);grid;hold on;(3)Notch 嵌波滤波器clear;Gb2=0.5;w0=0.35*pi;deltaw=0.1*pi;b=1/(1+tan(deltaw/2)*(sqrt(1-Gb2)/sqrt(Gb2))); B=[1 -2*cos(w0) 1].*b;A=[1 -2*b*cos(w0) (2*b-1)];w=0:pi/500:pi;H=freqz(B,A,w);plot(w/pi,abs(H));grid;(4)Peak 滤波器clear;Ac=3;Gb2=10^(-Ac/10);w0=0.35*pi;deltaw=0.1*pi;b=1/(1+tan(deltaw/2)*(sqrt(Gb2)/sqrt(1-Gb2))); B=[1 0 -1].*(1-b);A=[1 -2*b*cos(w0) (2*b-1)];w=0:pi/500:pi;H=freqz(B,A,w);plot(w/pi,abs(H));grid;(5)IIR低通滤波(Butterworth)% IIR Lowpass Use Butterworthclear;fs=20;fpass=4;fstop=5;Ap=0.5;As=10;wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;omegap=tan(wp/2);omegas=tan(ws/2);ep=sqrt(10^(Ap/10)-1);es=sqrt(10^(As/10)-1);N=ceil(log(es/ep)/log(omegas/omegap));omega0=omegap/ep^(1/N);K=floor(N/2);for i=1:Ktheta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:KG(i)=omega0^2/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka1(i)=2*(omega0^2-1)/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka2(i)=(1+2*omega0*cos(theta(i))+omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2); endif K<(N/2)G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);endw=0:pi/300:pi;Hw2=1./(1+(tan(w/2)/omega0).^(2*N));plot(w/pi,Hw2);grid;(6)IIR高通滤波(Butterworth)% IIR Hightpass Use Butterworthclear;fs=20;fpass=5;fstop=4;Ap=0.5;As=10;wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;omegap=cot(wp/2);omegas=cot(ws/2);ep=sqrt(10^(Ap/10)-1);es=sqrt(10^(As/10)-1);N=ceil(log(es/ep)/log(omegas/omegap));omega0=omegap/ep^(1/N);K=floor(N/2);for i=1:Ktheta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:KG(i)=omega0^2/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka1(i)=-2*(omega0^2-1)/(1-2*omega0*cos(theta(i))+omega0^2); endfor i=1:Ka2(i)=(1+2*omega0*cos(theta(i))+omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2); endif K<(N/2)G0=omega0/(omega0+1);a0=-(omega0-1)/(omega0+1);endw=(0+eps):pi/300:pi;Hw2=1./(1+(cot(w/2)/omega0).^(2*N));plot(w/pi,Hw2);grid;(7)IIR带通滤波(Butterworth)% IIR Bandpass Use Butterworthclear;fs=20;fpa=2;fpb=4;fsa=1.5;fsb=4.5;Ap=0.0877;As=16.9897;wpa=2*pi*fpa/fs;wpb=2*pi*fpb/fs;wsa=2*pi*fsa/fs;wsb=2*pi*fsb/fs;c=sin(wpa+wpb)/(sin(wpa)+sin(wpb));omegap=abs((c-cos(wpb))/sin(wpb));omegasa=(c-cos(wsa))/sin(wsa);omegasb=(c-cos(wsb))/sin(wsb);omegas=min(abs(omegasa),abs(omegasb));ep=sqrt(10^(Ap/10)-1);es=sqrt(10^(As/10)-1);N=ceil(log(es/ep)/log(omegas/omegap));omega0=omegap/ep^(1/N);K=floor(N/2);for i=1:Ktheta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:KG(i)=omega0^2/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka1(i)=4*c*(omega0*cos(theta(i))-1)/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka2(i)=2*(2*c^2+1-omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka3(i)=-(4*c*(omega0*cos(theta(i))+1))/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka4(i)=(1+2*omega0*cos(theta(i))+omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2); endG0=omega0/(1+omega0);a0(1)=-2*c/(1+omega0);a0(2)=(1-omega0)/(1+omega0); endw=(0+eps):pi/300:pi;Hw2=1./(1+((c-cos(w))./(omega0*sin(w))).^(2*N));plot(w/pi,Hw2);grid;(8)IIR带阻滤波(Butterworth)% IIR Bandstop Use Butterworthclear;fs=20;fpa=1.5;fpb=4.5;fsa=2;fsb=4;Ap=0.5;As=10;wpa=2*pi*fpa/fs;wpb=2*pi*fpb/fs;wsa=2*pi*fsa/fs;wsb=2*pi*fsb/fs;c=sin(wpa+wpb)/(sin(wpa)+sin(wpb));omegap=abs(sin(wpb)/(c-cos(wpb)));omegasa=sin(wsa)/(cos(wsa)-c);omegasb=sin(wsb)/(cos(wsb)-c);omegas=min(abs(omegasa),abs(omegasb));ep=sqrt(10^(Ap/10)-1);es=sqrt(10^(As/10)-1);N=ceil(log(es/ep)/log(omegas/omegap));omega0=omegap/ep^(1/N);K=floor(N/2);theta=zeros(1,K);theta(i)=pi*(N-1+2*i)/(2*N);endG=zeros(1,K);a1=zeros(1,K);a2=zeros(1,K);for i=1:KG(i)=omega0^2/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka1(i)=2*(omega0^2-1)/(1-2*omega0*cos(theta(i))+omega0^2);endfor i=1:Ka2(i)=(1+2*omega0*cos(theta(i))+omega0^2)/(1-2*omega0*cos(theta(i))+omega0^2); endif K<(N/2)G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);endw=(0+eps):pi/300:pi;Hw2=1./(1+(sin(w)./(omega0*(c-cos(w)))).^(2*N));plot(w/pi,Hw2);grid;(9)IIR低通滤波(chebyshev 1)% IIR Lowpass Use Chebyshev Type 1clear;fs=20;fpass=4;fstop=5;Ap=0.5;As=10;wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;omegap=tan(wp/2);omegas=tan(ws/2);ep=sqrt(10^(Ap/10)-1);es=sqrt(10^(As/10)-1);e=es/ep;w=omegas/omegap;N=ceil(log(e+sqrt(e^2-1))/log(w+sqrt(w^2-1)));a=log(1/ep+sqrt(1/ep^2+1))/N;omega0=omegap*sinh(a);K=floor(N/2);theta=zeros(1,K);omega=zeros(1,K);for i=1:Ktheta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:Komega(i)=omegap*sin(theta(i));endG=zeros(1,K);a1=zeros(1,K);a2=zeros(1,K);for i=1:KG(i)=(omega0^2+omega(i)^2)/(1-2*omega0*cos(theta(i))+omega0^2+omega(i)^2);endfor i=1:Ka1(i)=2*(omega0^2+omega(i)^2-1)/(1-2*omega0*cos(theta(i))+omega0^2+omega(i)^2) ;endfor i=1:Ka2(i)=(1+2*omega0*cos(theta(i))+omega0^2+omega(i)^2)/(1-2*omega0*cos(theta(i))+ omega0^2+omega(i)^2);endif K<(N/2)G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);elseH0=sqrt(1/(1+ep^2));endf=0:1/300:10;Hf2=1./(1+ep^2*(cheby(N,tan(pi*f/fs)/omegap)).^2);plot(f,abs(Hf2));grid;(9)IIR低通滤波(chebyshev 1)% IIR Lowpass Use Chebyshev Type 2clear;fs=20;fpass=4;fstop=5;Ap=0.5;As=10;wp=2*pi*fpass/fs;ws=2*pi*fstop/fs;omegap=tan(wp/2);omegas=tan(ws/2);ep=sqrt(10^(Ap/10)-1);es=sqrt(10^(As/10)-1);e=es/ep;w=omegas/omegap;N=ceil(log(e+sqrt(e^2-1))/log(w+sqrt(w^2-1)));a=log(es+sqrt(es^2+1))/N;omega0=omegas/sinh(a);K=floor(N/2);for i=1:Ktheta(i)=pi*(N-1+2*i)/(2*N);endfor i=1:Komega(i)=omegas/sin(theta(i));endfor i=1:KG(i)=(1+omega(i)^-2)/(1-2*omega0^-1*cos(theta(i))+omega0^-2+omega(i)^-2);endfor i=1:Ka1(i)=2*(1-omega0^-2+omega(i)^-2)/(1-2*omega0^-1*cos(theta(i))+omega0^-2+omeg a(i)^-2);endfor i=1:Ka2(i)=(1+2*omega0^-1*cos(theta(i))+omega0^-2+omega(i)^-2)/(1-2*omega0^-1*cos(th eta(i))+omega0^-2+omega(i)^-2);endfor i=1:Kb1(i)=2*(1-omega(i))/(1+omega(i));endif K<(N/2)G0=omega0/(omega0+1);a0=(omega0-1)/(omega0+1);elseH0=sqrt(1/(1+ep^2));endf=(0+eps):1/100:10;Hf2=(cheby(N,omegas./tan(pi*f/fs))).^2./((cheby(N,omegas./tan(pi*f/fs))).^2+es^2);plot(f,abs(Hf2));grid;(10)chebyshev 中用到的函数cheby.mfunction CN=cheby(N,x)if x<=1CN=cos(N*acos(x));elseCN=cosh(N*log(x+sqrt(x.^2-1)));end。
matlab中bode函数源程序 -回复

matlab中bode函数源程序-回复什么是Bode函数?在信号处理和控制系统中,Bode函数是一种常见的频率响应分析方法。
它用来描述线性时不变系统对输入信号频率的响应。
Bode函数由两个部分组成:振幅响应和相位响应。
振幅响应描述了系统对不同频率输入信号的增益,而相位响应描述了系统对输入信号的相位延迟。
Bode函数是通过将系统的传输函数表示为频率响应函数的形式而得到的。
传输函数是一个复数函数,它的实部表示系统的振幅响应,而虚部表示相位响应。
Bode函数将传输函数的振幅和相位部分分别转换为以频率为自变量的函数,从而更直观地观察系统对不同频率信号的响应。
在MATLAB中,Bode函数是一个内置函数,用于计算给定传输函数的振幅和相位响应。
它可以使用以下语法进行调用:[bode_mag, bode_phase, wout] = bode(num, den, w)其中,num是传输函数的分子多项式系数,den是传输函数的分母多项式系数,w是一个包含频率值的向量。
bode_mag、bode_phase和wout 分别是计算得到的振幅响应、相位响应和对应的频率输出。
使用Bode函数进行频率响应分析的步骤如下:1. 确定系统的传输函数:传输函数是一个将输入信号映射到输出信号的函数。
它由线性时不变系统的差分方程或微分方程表示。
传输函数可以用有理多项式表示,其中分子多项式表示输入信号的影响,分母多项式表示输出信号的影响。
2. 将传输函数转化为频率响应函数:使用MATLAB的tf函数将传输函数表示为频率域中的函数。
传输函数的分子多项式系数和分母多项式系数作为tf函数的输入。
3. 定义频率范围:确定频率范围,以便观察系统对不同频率的响应。
可以使用logspace函数生成一系列对数分布的频率值。
4. 使用Bode函数计算振幅和相位响应:使用Bode函数来计算给定传输函数的振幅和相位响应。
传输函数的分子多项式系数、分母多项式系数和频率向量作为输入。
数学建模2020c题matlab源程序

数学建模2020c题matlab源程序【实用版】目录一、数学建模 2020c 题概述二、Matlab 在数学建模中的应用三、2020c 题 matlab 源程序解析四、结论正文一、数学建模 2020c 题概述数学建模是一种运用数学方法和技术来解决实际问题的过程,其目的是通过建立数学模型,揭示问题的本质和规律,从而为解决实际问题提供理论依据。
2020c 题是数学建模竞赛中的一道题目,它涉及到了复杂的数学方法和技术,需要参赛者运用深厚的数学功底和灵活的思维方式来完成。
二、Matlab 在数学建模中的应用Matlab 是一种功能强大的数学软件,它既可以用于数学计算,也可以用于数据分析和可视化。
在数学建模过程中,Matlab 可以提供丰富的工具箱和函数,帮助用户解决各种复杂的数学问题。
例如,在 2020c 题中,Matlab 可以用于构建和求解数学模型,进行数据分析和绘图,从而为解题提供强大的支持。
三、2020c 题 matlab 源程序解析2020c 题的 matlab 源程序主要包括以下几个部分:(1)问题分析:这部分主要通过对题目的阅读和理解,明确问题的背景和要求,为构建数学模型打下基础。
(2)模型构建:这部分是数学建模的核心,需要根据问题的实际情况,构建出一个符合题意的数学模型。
在 2020c 题中,我们需要建立一个描述人口增长和资源消耗关系的模型。
(3)模型求解:这部分主要是运用 Matlab 的数学工具箱和函数,对建立的数学模型进行求解。
在 2020c 题中,我们需要运用 Matlab 的符号运算和数值计算功能,求解模型中的参数,从而得到人口增长和资源消耗的关系。
(4)结果分析:这部分主要是对模型求解的结果进行分析,得出问题的解决方案。
在 2020c 题中,我们需要根据模型求解的结果,分析人口增长和资源消耗的关系,从而为解决实际问题提供理论依据。
四、结论总的来说,数学建模 2020c 题是一道需要运用深厚数学功底和灵活思维方式的题目。
(完整版)遗传算法matlab实现源程序

附页:一.遗传算法源程序:clc; clear;population;%评价目标函数值for uim=1:popsizevector=population(uim,:);obj(uim)=hanshu(hromlength,vector,phen);end%obj%min(obj)clear uim;objmin=min(obj);for sequ=1:popsizeif obj(sequ)==objminopti=population(sequ,:);endendclear sequ;fmax=22000;%==for gen=1:maxgen%选择操作%将求最小值的函数转化为适应度函数for indivi=1:popsizeobj1(indivi)=1/obj(indivi);endclear indivi;%适应度函数累加总合total=0;for indivi=1:popsizetotal=total+obj1(indivi);endclear indivi;%每条染色体被选中的几率for indivi=1:popsizefitness1(indivi)=obj1(indivi)/total;endclear indivi;%各条染色体被选中的范围for indivi=1:popsizefitness(indivi)=0;for j=1:indivifitness(indivi)=fitness(indivi)+fitness1(j);endendclear j;fitness;%选择适应度高的个体for ranseti=1:popsizeran=rand;while (ran>1||ran<0)ran=rand;endran;if ran〈=fitness(1)newpopulation(ranseti,:)=population(1,:);elsefor fet=2:popsizeif (ran〉fitness(fet—1))&&(ran<=fitness(fet))newpopulation(ranseti,:)=population(fet,:);endendendendclear ran;newpopulation;%交叉for int=1:2:popsize-1popmoth=newpopulation(int,:);popfath=newpopulation(int+1,:);popcross(int,:)=popmoth;popcross(int+1,:)=popfath;randnum=rand;if(randnum〈 P>cpoint1=round(rand*hromlength);cpoint2=round(rand*hromlength);while (cpoint2==cpoint1)cpoint2=round(rand*hromlength);endif cpoint1>cpoint2tem=cpoint1;cpoint1=cpoint2;cpoint2=tem;endcpoint1;cpoint2;for term=cpoint1+1:cpoint2for ss=1:hromlengthif popcross(int,ss)==popfath(term)tem1=popcross(int,ss);popcross(int,ss)=popcross(int,term);popcross(int,term)=tem1;endendclear tem1;endfor term=cpoint1+1:cpoint2for ss=1:hromlengthif popcross(int+1,ss)==popmoth(term)tem1=popcross(int+1,ss);popcross(int+1,ss)=popcross(int+1,term);popcross(int+1,term)=tem1;endendclear tem1;endendclear term;endclear randnum;popcross;%变异操作newpop=popcross;for int=1:popsizerandnum=rand;if randnumcpoint12=round(rand*hromlength);cpoint22=round(rand*hromlength);if (cpoint12==0)cpoint12=1;endif (cpoint22==0)cpoint22=1;endwhile (cpoint22==cpoint12)cpoint22=round(rand*hromlength);if cpoint22==0;cpoint22=1;endendtemp=newpop(int,cpoint12);newpop(int,cpoint12)=newpop(int,cpoint22);newpop(int,cpoint22)=temp;。
919324-MATLAB程序设计与应用(第3版)-源程序-第12章 MATLAB Simulink系统仿真__源程序

第12章MATLAB Simulink系统仿真例12-8采用S函数实现y = nx,即把一个输入信号放大n倍。
①利用MATLAB语言编写S函数,程序如下。
%******************************************************* %S函数timesn.m,其输出是输入的n倍%******************************************************* function [sys,x0,str,ts]=timesn(t,x,u,flag,n)switch flagcase 0[sys,x0,str,ts]=mdlInitializeSizes; %初始化case 3sys=mdlOutputs(t,x,u,n); %计算输出量case {1,2,4,9}sys=[];otherwise %出错处理error(num2str(flag))end%******************************************************* %mdlInitializeSizes:当flag为0时进行整个系统的初始化%******************************************************* function [sys,x0,str,ts]=mdlInitializeSizes()%调用函数simsizes以创建结构sizessizes=simsizes;%用初始化信息填充结构sizessizes.NumContStates=0; %无连续状态sizes.NumDiscStates=0; %无离散状态sizes.NumOutputs=1; %有一个输出量sizes.NumInputs=1; %有一个输入信号sizes.DirFeedthrough=1; %输出量中含有输入量sizes.NumSampleTimes=1; %单个采样周期%根据上面的设置设定系统初始化参数sys=simsizes(sizes);%给其他返回参数赋值x0=[]; %设置初始状态为零状态str=[]; %将str变量设置为空字符串ts=[-1,0]; %假定继承输入信号的采样周期%初始化子程序结束%******************************************************* %mdlOutputs:当flag值为3时,计算输出量1%******************************************************* function sys=mdlOutputs(t,x,u,n)sys=n*u;%输出量计算子程序结束例12-9采用S函数来构造非线性分段函数。
Matlab源程序代码

正弦波的源程序:(一),用到的函数1,f2t函数function x=f2t(X)global dt df t f T N%x=f2t(X)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同并为2的整幂%本函数需要一个全局变量dt(时域取样间隔) X=[X(N/2+1:N),X(1:N/2)];x=ifft(X)/dt;end2,t2f函数。
function X=t2f(x)global dt df N t f T%X=t2f(x)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同,并为2的整幂。
%本函数需要一个全局变量dt(时域取样间隔) H=fft(x);X=[H(N/2+1:N),H(1:N/2)]*dt;end(二),主程序。
1,%(1)绘出正弦信号波形及频谱global dt df t f Nclose allk=input('取样点数=2^k, k取10摆布');if isempty(k), k=10; endf0=input('f0=取1(kz)摆布');if isempty(f0), f0=1; endN=2^k;dt=0.01; %msdf=1/(N*dt); %KHzT=N*dt; %截短期Bs=N*df/2; %系统带宽f=[-Bs+df/2:df:Bs]; %频域横坐标t=[-T/2+dt/2:dt:T/2]; %时域横坐标s=sin(2*pi*f0*t); %输入的正弦信号S=t2f(s); %S是s的傅氏变换a=f2t(S); %a是S的傅氏反变换a=real(a);as=abs(S);subplot(2,1,1) %输出的频谱plot(f,as,'b');gridaxis([-2*f0,+2*f0,min(as),max(as)])xlabel('f (KHz)')ylabel('|S(f)| (V/KHz)') %figure(2)subplot(2,1,2)plot(t,a,'black') %输出信号波形画图gridaxis([-2/f0,+2/f0,-1.5,1.5])xlabel('t(ms)')ylabel('a(t)(V)')gtext('频谱图')最佳基带系统的源程序:(一),用到的函数f2t函数和t2f函数。
(完整版)潮流计算的MATLAB源程序

%简单潮流计算的牛顿拉夫逊程序,相关的原始数据数据数据输入格式如下:%B1是支路参数矩阵,第一列和第二列是节点编号。
节点编号由小到大编写%对于含有变压器的支路,第一列为低压侧节点编号,第二列为高压侧节点%编号,将变压器的串联阻抗置于低压侧处理。
%第三列为支路的串列阻抗参数。
%第四列为支路的对地导纳参数。
%第五列为含变压器支路的变压器的变比%第六列为变压器是否含有变压器的参数,其中“1"为含有变压器,%“0”为不含有变压器。
%B2为节点参数矩阵,其中第一列为节点注入发电功率参数;第二列为节点负荷功率参数;第三列为节点电压参数;第六列为节点类型参数,其中“1”为平衡节点,“2”为PQ节点,“3”为PV节点参数。
%X为节点号和对地参数矩阵.其中第一列为节点编号,第二列为节点对地%参数。
n=input(’请输入节点数:n=');n1=input(’请输入支路数:n1=');isb=input(’请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input(’请输入支路参数:B1=’);B2=input(’请输入节点参数:B2=’);X=input('节点号和对地参数:X=');Y=zeros(n);Times=1; %置迭代次数为初始值%创建节点导纳矩阵for i=1:n1if B1(i,6)==0 %不含变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0。
5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0。
5*B1(i,4);else %含有变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)—1/(B1(i,3)*B1(i,5));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));endendYOrgS=zeros(2*n—2,1);DetaS=zeros(2*n—2,1); %将OrgS、DetaS初始化%创建OrgS,用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h—1,1)+real(B2(i,3))*(real(Y (i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))—imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:n %对PV节点的处理,注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h—1,1)=OrgS(2*h—1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))—real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendOrgS%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h—1,1);t=0;for i=1:nif B2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endendPVU%创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2h=h+1;DetaS(2*h—1,1)=real(B2(i,2))-OrgS(2*h—1,1);DetaS(2*h,1)=imag(B2(i,2))—OrgS(2*h,1);endendt=0;for i=1:n %对PV节点的处理,注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2—imag(B2(i,3))^2;endendDetaS%创建I,用于存储节点电流参数i=zeros(n—1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h—1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endendI%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n—2);h=0;k=0;for i=1:n %对PQ节点的处理if B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k—1)=—imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real (I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)—2*imag (I(h,1));else %非对角元素的处理Jacbi(2*h—1,2*k-1)=—imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h—1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=—Jacbi(2*h—1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h—1,2*k—1);endif k==(n-1) %将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0;endendendendendk=0;for i=1:n %对PV节点的处理if B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h—1,2*k-1)=—imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h—1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非对角元素的处理Jacbi(2*h-1,2*k—1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h—1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k—1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环的指针置于初始值,以确保雅可比矩阵换行k=0;endendendendendJacbi%求解修正方程,获取节点电压的不平衡量DetaU=zeros(2*n—2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理if B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j—1,1)*sqrt(—1);endendfor i=1:n %对PV节点的处理if B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB2%开始循环**********************************************************************while abs(max(DetaU))>prOrgS=zeros(2*n-2,1); %初始功率参数在迭代过程中是不累加的,所以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h—1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))—imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))—imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:nif i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real (B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))—real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendOrgS%创建DetaSh=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))—OrgS(2*h—1,1);DetaS(2*h,1)=imag(B2(i,2))—OrgS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h—1,1)=real(B2(i,2))—OrgS(2*h—1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;endendDetaS%创建Ii=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h—1,1)-OrgS(2*h,1)*sqrt(—1))/conj(B2(i,3));endendI%创建JacbiJacbi=zeros(2*n—2);h=0;k=0;for i=1:nif B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h—1,2*k—1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k—1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k—1)—2*imag(I(h,1));elseJacbi(2*h—1,2*k—1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h—1,2*k)=real(Y(i,j))*real (B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k—1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h—1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h—1,2*k-1)=—imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));elseJacbi(2*h-1,2*k—1)=-imag(Y(i,j))*real (B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n—1)k=0;endendendendendJacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:nif B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j—1,1)*sqrt(-1);endendfor i=1:nif B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(—1);endendB2Times=Times+1; %迭代次数加1endTimes一个原始数据的例子节点数 5支路数 5平衡节点编号 5精度pr 0。
lagrange插值matlab源程序

lagrange插值matlab源程序function y = lagrange_interpolation(x, y, x0)% 检查输入值的维度是否一致if numel(x) ~= numel(y)error('输入的x和y的长度不一致');end% 初始化插值多项式的系数n = length(x);p = zeros(n, n);p(:, 1) = 1; % p0 = 1% 计算拉格朗日插值基函数for i = 1:nfor j = 1:nif i ~= jp(i, j) = (x0 - x(j)) / (x(i) - x(j));endendend% 计算插值结果y = p * y; % 使用矩阵乘法计算插值结果end在这个函数中,输入参数x和y是已知的数据点,x0是需要计算插值的点。
然后,我们使用Lagrange插值基函数来计算插值多项式,并将结果存储在变量y中。
最后,返回y作为插值结果。
% 已知数据点x = [0, 1, 2, 3];y = [1, 2, 1, 3];% 需要计算插值的点x0 = 1.5;% 使用Lagrange插值进行估算y0 = lagrange_interpolation(x, y, x0);disp(y0); % 输出插值结果在这个例子中,我们使用了一组已知的数据点(x, y),然后使用Lagrange 插值方法来估算在x0 = 1.5处的y值。
然后,我们打印出这个估算结果。
请注意,这是一个简单的例子,实际使用时需要根据具体的数据和需求进行调整。
matlab中bode函数源程序

一、概述在控制系统工程中,频率响应是系统性能分析的重要手段之一。
Bode 图是频率响应的常用图示方法之一,它能够直观地展现系统的幅频特性和相频特性。
在MATLAB中,我们可以利用bode函数来绘制系统的Bode图,对系统的频率响应进行分析和评估。
二、bode函数的基本语法MATLAB中bode函数的基本语法如下:[bode_mag, bode_phase, w] = bode(sys)其中,sys表示系统的传递函数模型或状态空间模型;bode_mag和bode_phase分别表示系统的幅频特性和相频特性;w表示频率范围。
三、bode函数的使用方法1. 导入系统模型在使用bode函数之前,首先需要导入系统的传递函数模型或状态空间模型。
对于传递函数模型G(s),可以使用以下命令进行导入:sys = tf([1],[1 2 1])2. 绘制Bode图一旦导入了系统模型,就可以利用bode函数来绘制系统的Bode图。
使用以下命令可以实现:[bode_mag, bode_phase, w] = bode(sys)3. 显示Bode图绘制Bode图之后,可以使用以下命令来显示幅频特性和相频特性:figuresubplot(2,1,1)semilogx(w,20*log10(bode_mag))grid onxlabel('Frequency (rad/s)')ylabel('Magnitude (dB)')title('Bode Magnitude Plot')subplot(2,1,2)semilogx(w,bode_phase)grid onxlabel('Frequency (rad/s)')ylabel('Phase (deg)')title('Bode Phase Plot')四、实例演示下面我们以一个具体的系统为例,演示bode函数的使用方法。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
if (length(find(P<=0))~=0)
error('Not a prob.vector,negative component'); %判断是否符合概率分布条件
Huffman编码(1)
% huffman编码生成器%
%函数说明:%
% [W,L,q]=huffman(P)为huffman编码函数%
% P为信源的概率矢量,W为编码返回的码字%
% L为编码返回的平均码字长度,q为编码效率%
%*****************************************%
end
if (abs(sum(p)-1)>10e-10)
error('Not a prob.vector,component do not add up to 1') %判断是否符合概率和为1
end
% 1)排序
n=length(p);
x=1:n;
[p,x]=array(p,x);
% 2)计算代码组长度l
end
% 1)初始化Pa
[r,s]=size(P);
Pa=(1/(r+eps))*ones(1,r);
sumrow=zeros(1,r);
Pba=P;
% 2)进行迭代计算
n=0;
C=0;
CC=1;
while abs(CC-C)>=k
n=n+1;
% (1)先求Pb
Pb=zeros(1,s);
for j=1:s
for i=1:r
Pb(j)=Pb(j)+Pa(i)*Pba(i,j);
end
end
% (2)再求Pab
suma=zeros(1,s);
for j=1:s
for i=1:r
Pab(j,i)=Pa(i)*Pba(i,j)/(Pb(j)+eps);
suma(j)=suma(j)+Pa(i)*Pba(i,j)*log2((Pab(j,i)+eps)/(Pa(i)+eps));
end
% 3)求信道容量C
C=sum(suma);
% 4)求下一次Pa,即Paa
L=zeros(1,r);
sumaa=0;
for i=1:r
for j=1:s
L(i)=L(i)+Pba(i,j)*log(Pab(j,i)+eps);
end
a(i)=exp( L(i));
end
sumaa=sum(a);
n=length(current_P);
add(1)=current_P(1);
% add(i)=0;
add(i)=add(i-1)+current_P(i);
end
% 2)求概率和最接近的两小组
s=add(n);
for i=1:n
temp(i)=abs(s-2*add(i));
% next_index为返回的下次分组的信源的下标%
% code_num为返回的ASCII值%
%*********************************************************************%
function [next_P,code_num,next_index]=compare(current_P,current_index);
if (i<n)
for k=(maxN-1):-1:i
P(:,k+1)=P(:,k);
end
end
end
P(:,i)=MAX;
end
p=P(1,:);
x=P(2,:);
% shannon编码生成器%
%函数说明:%
% [W,L,q]=shannon(p)为shannon编码函数%
% p为信源的概率矢量,W为编码返回的码字%
function [W,L,q]=fano(P)
%提示错误信息
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') %判断是否符合概率和为1
end
% 1)排序
n=length(P);
x=1:n;
[P,x]=array(P,x);
% 2)将信源符号分组并得到对应的码字
for i=1:n
current_index=i;
j=1;
current_P=P;
while 1
[next_P,code_num,next_index]=compare(current_P,current_index);
% Pb:输出概率矩阵,Pab:反向转移概率矩阵%
% C:初始信道容量,r:输入符号数,s:输出符号数%
%**************************************************%
function [CC,Paa]=ChannelCap(P,k)
%提示错误信息
if (length(find(P<0))~=0)
next_P=current_P((k+1):n);
end
% fano编码生成器%
%函数说明:%
% [W,L,q]=fano(P)为fano编码函数%
% P为信源的概率矢量,W为编码返回的码字%
% L为编码返回的平均码字长度,q为编码效率%
%*****************************************%
current_index=next_index;
current_P=next_P;
W(i,j)=code_num;
j=j+1;
if (length(current_P)==1)
break;
end
end
l(i)=length(find(abs(W(i,:))~=0)); %得到各码字的长度
end
L=sum(P.*l); %计算平均码字长度
P(i)=P(i)*2-temp(i,j);
end
end
% b)给W赋ASCII码值,用于显示二进制代码组W
for i=1:n
for j=1:l(i)
if (temp(i,j)==0)
W(i,j)=48;
else
W(i,j)=49;
end
end
end
L=sum(p.*l); %计算平均码字长度
H=entropy(p,2); %计算信源熵
附录B 离散无记忆信道容量的迭代计算
%信道容量C的迭代算法%
%函数说明:%
% [CC,Paa]=ChannelCap(P,k)为信道容量函数%
%变量说明:%
% P:输入的正向转移概率矩阵,k:迭代计算精度%
% CC:最佳信道容量,Paa:最佳输入概率矩阵%
% Pa:初始输入概率矩阵,Pba:正向转移概率矩阵%
end
[c,k]=min(temp);
% 3)对分组的信源赋ASCII值
if (current_index<=k)
next_index=current_index;
code_num=48;
next_P=current_P(1:k);
else
next_index=current_index-k;
code_num=49;
% L为编码返回的平均码字长度,q为编码效率%
%*****************************************%
function [W,L,q]=shannon(p)
%提示错误信息
if (length(find(p<=0))~=0)
error('Not a prob.vector,negative component'); %判断是否符合概率分布条件
s0='很好!输入正确,编码结果如下:';
s1='Fano编码所得码字W:';
s2='Fano编码平均码字长度L:';
s3='Fano编码的编码效率q:';
disp(s0);
disp(s1),disp(B),disp(W);
disp(s2),disp(L);
disp(s3),disp(q);
附录E Huffman编码
% [next_P,next_index,code_num]=compare(current_P,current_index)%
%为比较函数,主要用于信源符号的分组%
% current_P为当前分组的信源的概率矢量%
% current_index为当前分组的信源的下标%
% next_P为返回的下次分组的信源的概率矢量%
disp(s3),disp(n);
附录C Shannon编码
%函数说明:%
% [p,x]=array(P,X)为按降序排序的函数%
% P为信源的概率矢量,X为概率元素的下标矢量%
% p为排序后返回的信源的概率矢量%
% x为排序后返回的概率元素的下标矢量%
%*******************************************%
H=entropy(P,2); %计算信源熵