Matlab源程序代码
复杂网络聚类系数和平均路径长度计算的MATLAB源代码
复杂网络聚类系数和平均路径长度计算的MATLAB源代码复杂网络的聚类系数和平均路径长度是度量网络结构特征的重要指标。
下面是MATLAB源代码，用于计算复杂网络的聚类系数和平均路径长度。
首先，我们需要定义一个函数，用于计算节点的聚集系数。
这个函数的输入参数是邻接矩阵和节点的索引，输出参数是节点的聚类系数。
```matlabfunction cc = clustering_coefficient(adj_matrix, node_index) neighbors = find(adj_matrix(node_index, :));k = length(neighbors);if k < 2cc = 0;elseconnected_count = 0;for i = 1:k-1for j = i+1:kif adj_matrix(neighbors(i), neighbors(j))connected_count = connected_count + 1;endendendcc = 2 * connected_count / (k * (k - 1));endend```接下来，我们定义一个函数，用于计算整个网络的平均聚合系数。
```matlabfunction avg_cc = average_clustering_coefficient(adj_matrix) n = size(adj_matrix, 1);cc = zeros(n, 1);for i = 1:ncc(i) = clustering_coefficient(adj_matrix, i);endavg_cc = sum(cc) / n;end```然后，我们需要计算网络的平均最短路径长度。
这里我们使用了Floyd算法来计算每对节点之间的最短路径。
```matlabfunction avg_path_length =average_shortest_path_length(adj_matrix)n = size(adj_matrix, 1);dist_matrix =graphallshortestpaths(sparse(double(adj_matrix)));avg_path_length = sum(dist_matrix(:)) / (n^2 - n);end```最后，我们可以使用这些函数来计算一个复杂网络的聚类系数和平均路径长度。
matlab的conv的c源代码
matlab的conv的c源代码MATLAB的conv函数是非常常用的信号处理函数，它用于进行离散卷积运算。
在MATLAB中，conv函数的底层实现是使用C语言编写的，我们可以通过查看源代码来了解其具体实现细节。
以下是MATLAB的conv函数的部分C源代码：```c#include "mex.h"/* Gateway Function */void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){/* Check input and output arguments */if (nrhs != 2){mexErrMsgIdAndTxt("MATLAB:conv:invalidNumInputs","Two input arguments required.");}if (nlhs > 1){mexErrMsgIdAndTxt("MATLAB:conv:maxlhs","Too many output arguments.");}/* Get input data */double *input1 = mxGetPr(prhs[0]);double *input2 = mxGetPr(prhs[1]);mwSize len1 = mxGetNumberOfElements(prhs[0]);mwSize len2 = mxGetNumberOfElements(prhs[1]);/* Calculate output size */mwSize outlen = len1 + len2 - 1;/* Create output array */plhs[0] = mxCreateDoubleMatrix(1, outlen, mxREAL); double *output = mxGetPr(plhs[0]);/* Perform convolution */for (mwSize i = 0; i < outlen; i++){output[i] = 0;for (mwSize j = 0; j < len2; j++){if (i - j >= 0 && i - j < len1){output[i] += input1[i - j] * input2[j];}}}}```以上是MATLAB的conv函数的简化版本C源代码。
TDOA定位的Chan算法MATLAB源代码
TDOA定位的Chan算法MATLAB源代码func on [POS_ref,POS1,POS2,POS3,POS4] = TDOA_chan(R,Pbs,Q)%% TDOA定位定位Chan算法%*********************************************************% CHAN算法,假设移动台与各基站位置较近,需进行三四WLS计算% 输入参数:% R(N-1×1): TDOA测量值% Pbs(N×2): 基站的坐标,第一列为X,第二列为Y;参考基站坐标位于第一行% 输出参数:% POS_ref(2X1):第一次WLS的估计结果,作为参考作最终定位结果的判决% POS1(2X1) :定位结果1% POS2(2X1) :定位结果2% POS3(2X1) :定位结果3% POS4(2X1) :定位结果4%这段程序用于自己产生算法输入参数,用于算法测试% delta = 10; %TDOA测量误差的标准差,用于产生Q矩阵% M=4; %参考基站外的基站数量% Xb = 20;Yb = 100; %参考的基站坐标% X = zeros(M,2);% a = 2*pi/M;% for i=1:M %生成其他基站的坐标% X(i,1) = 400*cos(a*(i-1));% X(i,2) = 400*sin(a*(i-1));% end% Xreal = -150;Yreal = 200;%移动台真实坐标% %产生TDOA测量结果% Rb = sqrt((Xreal - Xb)^2+(Yreal - Yb)^2);%移动台到基站的真实距离% N = normrnd(0,delta,1,M);%产生TDOA测量误差% Kb = Xb^2+Yb^2;% R = zeros(M,1);% for i=1:M %产生TDOA测量值% R(i) = -Rb+sqrt((Xreal - X(i,1))^2+(Yreal - X(i,2))^2)+N(i);% end% Q = (0.5*eye(M)+0.5*ones(M))*(delta^2);% Pbs = [Xb Yb;X];% N = 5;% [POS_ref POS1 POS2 POS3 POS4] = TDOA_chan(R,Pbs,Q)N = size(Pbs,1);K = zeros(1,N);K = Pbs(:,1).^2 + Pbs(:,2).^2;ha = 0.5*(R.^2-K(2:N)+K(1))Ga = -[Pbs(2:N,1)-Pbs(1,1) Pbs(2:N,2)-Pbs(1,2) R]%计算第一次WLS估计结果（远距算法）Za = inv(Ga.'*inv(Q)*Ga)*Ga.'*inv(Q)*ha;Ba = diag(Va);Fa = Ba*Q*Ba;Zacov = inv(Ga.'*inv(Fa)*Ga);%第一次ＷＬＳ计算（近距算法）Va2 = sqrt((Za2(1)-Pbs(2:N,1)).^2+(Za2(2)-Pbs(2:N,2)).^2);Ba2 = diag(Va2);Fa2 = Ba2*Q*Ba2;Zacov2 = inv(Ga.'*inv(Fa2)*Ga);%第二次ＷＬＳ计算（近距算法）Gb = [1 0;0 1;1 1];Bb2 = [Za2(1)-Pbs(1,1) 0 0;0 Za2(2)-Pbs(1,2) 0;0 0 sqrt((Za2(1)-Pbs(1,1))^2+ (Za2(2)-Pbs(1,2))^2)];Fb2 = 4*Bb2*Zacov2*Bb2;POS_ref = Za2(1:2,1); %第一次WLS的估计结果,作为参考作最终定位结果的判决POS1 = sqrt(Zb2)+Pbs(1,:).'; %定位结果1POS2 = -sqrt(Zb2)+Pbs(1,:).';%定位结果2POS3 = [-sqrt(Zb2(1));sqrt(Zb2(2))]+Pbs(1,:).';%定位结果3POS4 = [sqrt(Zb2(1));-sqrt(Zb2(2))]+Pbs(1,:).';%定位结果4源代码运行结果展示037附图1.emf。
互信息matlab源代码
互信息matlab源代码互信息是一种用于评估两个变量之间相关性的统计量。
它可以用于特征选择、机器学习和数据挖掘等领域。
下面是一个用MATLAB编写的计算两个随机变量之间互信息的例子。
首先，我们定义两个随机变量X和Y。
这里我们使用正态分布生成随机数据。
```matlabmu1 = 0; mu2 = 0; sigma1 = 1; sigma2 = 1; rho = 0.5;n = 1000;接下来，我们定义一个函数`calculate_entropy`来计算随机变量的熵。
熵用于衡量随机变量的不确定性。
```matlabfunction entropy = calculate_entropy(data)frequency = tabulate(data); % 统计每个值的频率probability = frequency(:,3)/100; % 计算每个值的概率probability = probability(probability ~= 0); % 剔除概率为0的值entropy = -sum(probability.*log2(probability)); % 计算熵end```接下来，我们定义另一个函数`calculate_mutual_information`来计算X和Y之间的互信息。
```matlabfunction mi = calculate_mutual_information(x, y)joint_probability = zeros(length(unique(x)), length(unique(y)));for i=1:length(x)x_idx = find(unique(x) == x(i));y_idx = find(unique(y) == y(i));joint_probability(x_idx, y_idx) = joint_probability(x_idx, y_idx) + 1;endjoint_probability = joint_probability/sum(sum(joint_probability)); % 计算联合概率分布x_probability = sum(joint_probability,2); % 计算X的概率分布y_probability = sum(joint_probability,1); % 计算Y的概率分布entropy_x = calculate_entropy(x); % 计算X的熵entropy_y = calculate_entropy(y); % 计算Y的熵mi = 0;for i=1:size(joint_probability,1)for j=1:size(joint_probability,2)if joint_probability(i,j) ~= 0mi = mi +joint_probability(i,j)*log2(joint_probability(i,j)/(x_probability(i)*y_probabi lity(j)));endendendend```最后，我们使用这两个函数计算X和Y之间的互信息。
chen系统李雅普诺夫指数的matlab源代码
chen系统李雅普诺夫指数的matlab源代码以下是计算 Chen 系统的李雅普诺夫指数的 Matlab 源代码。
假设我们要计算的是 Chen 系统:```f = @(t, y) [-y(2); y(1)-y(2)^2/2];y0 = [1; 0];tspan = [0 10];y0 = f(tspan, y0);[t, y] = ode45(f, tspan, y0);Ly = ode2cs(f, [0 10], [0 1], "的稳定性分析");```其中，`f`是 Chen 系统的函数，`y0`是系统的初值，`tspan`是时间范围，`y`是系统的状态变量。
`Ly`是李雅普诺夫指数，它表示系统的稳定性。
Matlab 代码的具体步骤如下:1. 定义 Chen 系统的函数`f`,并设置初值`y0`。
2. 使用`ode45`函数求解 Chen 系统的 ODE。
3. 使用`ode2cs`函数计算李雅普诺夫指数`Ly`。
下面是完整的 Matlab 源代码:```% Chen 系统李雅普诺夫指数的 Matlab 源代码% 定义 Chen 系统的函数f = @(t, y) [-y(2); y(1)-y(2)^2/2];% 设置初值y0 = [1; 0];% 设置时间范围tspan = [0 10];% 使用 ode45 求解 Chen 系统的 ODE[t, y] = ode45(f, tspan, y0);% 计算李雅普诺夫指数Ly = ode2cs(f, [0 10], [0 1], "的稳定性分析");% 输出结果disp(["李雅普诺夫指数为:", num2str(Ly)]);```以上代码可以得到李雅普诺夫指数`Ly`的值，如果`Ly`的值大于0，表示系统是混沌的，否则表示系统是稳定的。
因子分析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代码示例，你可以根据自己的需求修改和扩展代码。
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 strcat函数源代码
MATLAB中的strcat函数是用来连接字符串的函数，它可以将多个字符串连接成一个新的字符串。
1. strcat函数的定义strcat函数的定义如下所示：```matlabfunction s = strcat(varargin)```其中，s为连接后的字符串，varargin为输入的字符串参数。
2. strcat函数的使用方法strcat函数的使用方法如下所示：```matlabs = strcat(str1, str2, ..., strN)```其中，str1、str2、...、strN为要连接的字符串参数，s为连接后的字符串。
3. strcat函数的源代码strcat函数的源代码如下所示：```matlabfunction s = strcat(varargin)narginchk(1, inf)for i = 1:narginvalidateattributes(varargin{i}, {'char'}, {'row'})ends = [varargin{:}];end```在上面的源代码中，首先使用了narginchk函数来检查输入参数的个数，确保参数个数在1到无穷大之间。
然后使用了validateattributes 函数来验证输入参数的类型，确保输入的参数都是字符型的向量。
使用了括号运算符[]来将所有的字符串参数连接成一个新的字符串，并将结果赋值给变量s。
4. strcat函数的示例下面是使用strcat函数的一个示例：```matlabstr1 = 'Hello, ';str2 = 'world!';s = strcat(str1, str2);disp(s);```执行上面的示例代码后，变量s的值为'Hello, world!'。
5. strcat函数的注意事项在使用strcat函数时，需要注意以下几点：- 输入的参数必须是字符型的向量，否则会导致错误。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
正弦波的源程序：(一)，用到的函数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函数。
代码>>(二),主程序globaldt t f df N Tclose allclear Eb_N0 Pek=input('取样点数＝2^k, k取13左右'); if isempty(k), k=13; endz=input('每个信号取样点数＝2^z, z<k'); if isempty(z), z=3; endaa=input('滚降系数=[0.5]');if aa==[],aa=0.5;endN=2^k;L=2^z;M=N/L;Rb=2; %码速率是2Mb/sTs=1/Rb; %码元间隔dt=Ts/L; %时域采样间隔df=1/(N*dt); %频域采样间隔T=N*dt; %截短时间Bs=N*df/2; %系统带宽f=[-Bs+df/2:df:Bs]; %频域横坐标t=[-T/2+dt/2:dt:T/2]; %时域横坐标g=sin(pi*t/Ts).*cos(pi*t*aa/Ts)./[pi*t/Ts.*(1-4*t.^2*aa^2/Ts^2)]; %升余弦脉冲波形GG=t2f(g);GG=abs(GG); %升余弦脉冲的傅式变换GT=sqrt(GG);GR=GT; %最佳系统的发送接收滤波器的傅式变换for l1=1:20;Eb_N0(l1)=(l1-1); %Eb/N0 in dBeb_n0(l1)=10^(Eb_N0(l1)/10);Eb=1;n0=Eb/eb_n0(l1); %信道的噪声谱密度sita=n0*Bs; %信道中噪声功率n_err=0; %误码计数for l2=1:20;b=sign(randn(1,M));s=zeros(1,N); %产生冲激序列s(L/2:L:N)=b/dt;SS=t2f(s);S=SS.*GG; %信道的傅式变化a=f2t(S);a=real(a); %不加噪声的输出n_ch=sqrt(sita)*randn(size(t)); %信道噪声N_CH=t2f(n_ch);nr=real(f2t(N_CH.*GR)); %输出噪声sr=a+nr; %接收信号y=sr(L/2-L/4:L:N-L/4); %取样bb=sign(y); %判决n_err=n_err+length(find(bb~=b));%错误累计endPe(l1)=n_err/(M*l2);semilogy(Eb_N0,Pe,'g'); %Pe~Eb/N0曲线画图xlabel('Eb/N0');ylabel('Pe');title('Pe~Eb/N0曲线');eb_n0=10.^(Eb_N0/10);hold onsemilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)));axis([0,15,1e-4,1])xlabel('Eb/N0')ylabel('Pe')legend('实际的','理论的')end多径传输的源程序：(一)，用到的函数f2t函数和t2f函数。
代码>>(二),主程序global dt t f df N Tclose allclear Eb_N0 Pek=input('取样点数＝2^k, k取13左右');if isempty(k), k=13; endz=input('每个信号取样点数＝2^z, z<k');if isempty(z), z=3; endaa=input('滚降系数=[0.5]');if aa==[],aa=0.5;endN=2^k;L=2^z;M=N/L;Rb=2; %码速率是2Mb/sTs=1/Rb; %码元间隔dt=Ts/L; %时域采样间隔df=1/(N*dt); %频域采样间隔T=N*dt; %截短时间Bs=N*df/2; %系统带宽f=[-Bs+df/2:df:Bs]; %频域横坐标t=[-T/2+dt/2:dt:T/2]; %时域横坐标g=sin(pi*t/Ts).*cos(pi*t*aa/Ts)./[pi*t/Ts.*(1-4*t.^2*aa^2/Ts^2)]; %升余弦脉冲波形GG=t2f(g);GG=abs(GG); %升余弦脉冲的傅式变换GT=sqrt(GG);GR=GT; %最佳系统的发送接收滤波器的傅式变换for l1=1:20;Eb_N0(l1)=(l1-1); %Eb/N0 in dBeb_n0(l1)=10^(Eb_N0(l1)/10);Eb=1;n0=Eb/eb_n0(l1); %信道的噪声谱密度sita=n0*Bs; %信道中噪声功率n_err=0; %误码计数for l2=1:5;dm=Ts/2;C=abs(1-0.5*exp(-j*2*pi*f*dm));%多径信道b=sign(randn(1,M));s=zeros(1,N); %产生冲激序列s(L/2:L:N)=b/dt;SS=t2f(s);S=SS.*GG.*C; %信道的傅式变化a=f2t(S);a=real(a); %不加噪声的输出n_ch=sqrt(sita)*randn(size(t)); %信道噪声N_CH=t2f(n_ch);nr=real(f2t(N_CH.*GR)); %输出噪声sr=a+nr; %接收信号y=sr(L/2:L:N); %取样bb=sign(y); %判决n_err=n_err+length(find(bb~=b));%错误累计endPe(l1)=n_err/(M*l2);semilogy(Eb_N0,Pe,'black'); %Pe~Eb/N0曲线画图xlabel('Eb/N0');ylabel('Pe');title('Pe~Eb/N0曲线'); eb_n0=10.^(Eb_N0/10);hold onsemilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)));axis([0,15,1e-3,1])xlabel('Eb/N0')ylabel('Pe')legend('实际的','理论的')end高斯噪声的源程序：function[gsrv1,gsrv2]=gngauss(m,sgma)if nargin==0,m=0;sgma=1;elseif nargin==1,sgma=m;m=0;end;u=rand;z=sgma*(sqrt(2*log(1/(1-u))));u=rand;gsrv1=m+z*cos(2*pi*u);gsrv2=m+z*sin(2*pi*u);Q函数互补的源程序：function[y]=Qfunct(x)y=(1/2)*erfc(x/sqrt(2));单极性归零码的源程序：(一)，用到的函数f2t函数和t2f函数。
代码>>(二),主程序global dt t f df N Tclose allk=input('取样点数＝2^k, k取10左右'); if isempty(k), k=13; endz=input('每个信号取样点数＝2^z,z<k'); if isempty(z), z=5; endf0=input('f0=取1(kz)左右');if isempty(f0), f0=1; endN=2^kL=2^z;M=N/L;Rb=2;Ts=0.5 %码元宽度是0.5usdt=Ts/L; %时域采样间隔df=1/(N*dt) %MHzT=N*dt %截短时间Bs=N*df/2 %系统带宽f=[-Bs+df/2:df:Bs]; %频域横坐标t=[-T/2+dt/2:dt:T/2]; %时域横坐标EP=zeros(size(f))+eps;for ii1=1:30b=round(rand(1,M)); %产生随机码s=zeros(1,N);for ii=1:L/2;s(ii+[0:M-1]*L)=b;endS=t2f(s); %S是s的傅氏变换a=f2t(S); %a是S的傅氏反变换a=abs(a);P=S.*conj(S)/T; %功率谱EP=(EP*(ii1-1)+P+eps)/ii1;figure(1) %输出的功率画图plot(f,10*log10(EP),'b');gridaxis([-15,15,-100,20]);xlabel('f (KHz)');ylabel('P (V/KHz)') ;title('功率谱图') figure(2) %输出信号波形画图plot(t,a,'b')gridaxis([-3,3,0,1.3]);xlabel('t (ms)');ylabel('s(t) (V)');title('单极性RZ的时域图') end8psk的源程序：function[pb,ps]=cm_sm32(snr_in_dB)N=10000;E=1;snr=10^(snr_in_dB/10);sgma=sqrt(E/(3*2*snr));s000=[1 0];s001=[cos(pi/4) sin(pi/4)];s011=[0 1];s010=[cos(3*pi/4) sin(3*pi/4)];s110=[-1 0];s111=[cos(5*pi/4) sin(5*pi/4)];s101=[0 -1];s100=[cos(7*pi/4) sin(7*pi/4)];for i=1:N,temp=rand;if(temp<0.125),dsource1(i)=0;dsource2(i)=0;dsource3(i)=0;elseif(temp<0.25),dsource1(i)=0;dsource2(i)=0;dsource3(i)=1;elseif(temp<0.375),dsource1(i)=0;dsource2(i)=1;dsource3(i)=0;elseif(temp<0.5),dsource1(i)=0;dsource2(i)=1;dsource3(i)=1;elseif(temp<0.625),dsource1(i)=1;dsource2(i)=0;dsource3(i)=0;elseif(temp<0.75),dsource1(i)=1;dsource2(i)=0;dsource3(i)=1;elseif(temp<0.875),dsource1(i)=1;dsource2(i)=1;dsource3(i)=0;elsedsource1(i)=1;dsource2(i)=1;dsource3(i)=1;end;end;numofsymbolerror=0;numofbiterror=0;for i=1:N,n(1)=gngauss(sgma);n(2)=gngauss(sgma);if((dsource1(i)==0)&(dsource2(i)==0)&(dsource3(i)==0)),r=s000+n;elseif((dsource1(i)==0)&(dsource2(i)==0)&(dsource3(i)==1)), r=s001+n;elseif((dsource1(i)==0)&(dsource2(i)==1)&(dsource3(i)==0)), r=s010+n;elseif((dsource1(i)==0)&(dsource2(i)==1)&(dsource3(i)==1)), r=s011+n;elseif((dsource1(i)==1)&(dsource2(i)==0)&(dsource3(i)==0)), r=s100+n;elseif((dsource1(i)==1)&(dsource2(i)==0)&(dsource3(i)==1)),r=s101+n;elseif((dsource1(i)==1)&(dsource2(i)==1)&(dsource3(i)==0)), r=s110+n;elser=s111+n;end;c000=dot(r,s000);c001=dot(r,s001);c010=dot(r,s010);c011=dot(r,s011);c100=dot(r,s100);c101=dot(r,s101);c110=dot(r,s110);c111=dot(r,s111);c_max=max([c000 c001 c010 c011 c100 c101 c110 c111]);if(c000==c_max),decis1=0;decis2=0;decis3=0;elseif(c001==c_max),decis1=0;decis2=0;decis3=1;elseif(c010==c_max),decis1=0;decis2=1;decis3=0;elseif(c011==c_max),decis1=0;decis2=1;decis3=1;elseif(c100==c_max),decis1=1;decis2=0;decis3=0;elseif(c101==c_max),decis1=1;decis2=0;decis3=1;elseif(c110==c_max),decis1=1;decis2=1;decis3=0;elsedecis1=1;decis2=1;decis3=1;end;symbolerror=0;if(decis1~=dsource1(i)),numofbiterror=numofbiterror+1;symbolerror=1;end;if(decis2~=dsource2(i)),numofbiterror=numofbiterror+1;symbolerror=1;end;if(decis3~=dsource3(i)),numofbiterror=numofbiterror+1;symbolerror=1;end;if(symbolerror==1),numofsymbolerror=numofsymbolerror+1;end;end;ps=numofsymbolerror/N;pb=numofbiterror/(3*N);4PSK的源程序：% 对M=4的PSK通信系统进行蒙特卡罗仿真%echo onSNRindB1=0:2:10;SNRindB2=0:0.1:10;for i=1:length (SNRindB1),[pb , ps]=cm_sm32 (SNRindB1 (i)); %仿真比特和符号误码率smld_bit_err_prb (i)=pb;smld_symbol_err_prb (i)=ps;end;for i=1:length (SNRindB2),SNR=exp(SNRindB2(i)*log(10)/10); %信噪比theo_err_prb(i)=Qfunct(sqrt(2*SNR)); %理论比特误码率end;%随后绘图曲线semilogy (SNRindB1,smld_bit_err_prb,'*');holdsemilogy (SNRindB1,smld_symbol_err_prb,'O');semilogy (SNRindB2,theo_err_prb);gridxlabel('Eb/No in dB')ylabel('error probability')title('4PSK通信系统的蒙特卡洛仿真')gtext('（注：“—”理论误码率；“*”误比特率；“o”误符号率）')。