MATLAB，3x3mm的光感接收器上，40X40的房间，每个房间（细胞）75X75，有效接收区域（活跃区域）65X65。已知光子落点坐标（x，y）有9000个。
三个要求
1号:不管活性区域,进入每个细胞的光子数。（落入75x75的光子数）
2号:进入每个细胞活性区域的光子数。（落入65X65的光子数）
3号:PD适用值(0.563)来获取并选择值
最后,为了做3号,要经过1号和2号。
下面的代码已经完成三号，我的问题是一号和二号要求代码怎么加入下面的代码里。我自己也写了一下不过我刚开始学，希望可以给我一个参考。
这里是示意图：
%%光子消减
clear all；clc；
％ 重置
数据 = cell2mat(readcell("COUNTS1"));
%读取数据
all_length = 3；gapd_length = 0.025；active_length = 0.015；inactive_length = 0.005；
% 单元格大小和活动区域设置
gapd_active_count = zeros([all_length/gapd_length all_length/gapd_length]);
%根据单元格大小创建数组
对于我 = 1：大小（数据，1）
% 重复进行到输PHOTON 的位置
temp_x = floor(data(i,1)/gapd_length);temp_y = floor(data(i,2)/gapd_length);
%用PHOTON 检测细胞位置
如果 ((data(i,1)-temp_x*gapd_length)>=inactive_length)&&((data(i,1)-temp_x*gapd_length)<=(inactive_length+active_length)) && ((data(i,2)- temp_y*gapd_length)>=inactive_length)&&((data(i,2)-temp_y*gapd_length)<=(inactive_length+active_length))
% CELL 判断是否为活动区域，进行条件语句
gapd_active_count(temp_y+1, temp_x+1) = gapd_active_count(temp_y+1, temp_x+1) + 0.563;
%如果为真，则添加 0.563 的 PD 值
end
gapd_active_count(gapd_active_count>1) = 1;
%由于值不能超过1，所以超过1就设置为1
总和（总和（gapd_active_count））
% 添加所有落入活动区域的值
这样应该可以吧
data = rand(10000,2)*3
[cell_count,area_count]=get_map(data)
h =heatmap(cell_count);
title(h,'cell count')
h = heatmap(area_count);
title(h,'area count')
function [cell_count,area_count]=get_map(data)
w = 0.075;
area_b = [0.005, w-0.005];
cell_count = zeros(40,40);
area_count = zeros(40,40);
re = mod(data,w);
location = int32((data-re)/w);
for i=1:size(location,1)
x = location(i,1)+1;
y = location(i,2)+1;
cell_count(x, y) = cell_count(x, y) + 1;
if all(re(i)>area_b(1)) && all(re(i)<area_b(2))
area_count(x,y) = area_count(x,y) + 1;
end
end
end
我会尝试给出一些可能的优化代码。 对于1号要求，我们只需要针对每个单元格计算光子的总数，而不需要考虑该光子是否在区域内。 这可以通过直接将所有光子分配到它们的单元格中并计算每个单元格的光子数量来完成。 对于2号要求，我们可以修改这种方法，仅计算位于活动区域内的光子数量。 这可以通过仅将位于活动区域内的光子分配到单元格中并计算每个单元格的光子数量来完成。 以下是可能的优化代码示例：
%%光子消减 clear all; clc; % 重置 data = cell2mat(readcell("COUNTS1")); %读取数据 all_length = 3; gapd_length = 0.025; active_length = 0.015; inactive_length = 0.005; % 单元格大小和活动区域设置 gapd_count = zeros([all_length/gapd_length all_length/gapd_length]); gapd_active_count = zeros([all_length/gapd_length all_length/gapd_length]); %根据单元格大小创建数组 for i = 1:size(data,1) % 分配光子到单元格中 temp_x = floor(data(i,1)/gapd_length);temp_y = floor(data(i,2)/gapd_length); if((temp_x > 0) && (temp_x < (all_length/gapd_length + 1)) && (temp_y > 0) && (temp_y < (all_length/gapd_length+1))) gapd_count(temp_x,temp_y) = gapd_count(temp_x,temp_y) + 1; if((data(i,1) - temp_xgapd_length - gapd_length/2) * (data(i,1) - temp_xgapd_length - gapd_length/2) + (data(i,2) - temp_ygapd_length - gapd_length/2) * (data(i,2) - temp_ygapd_length - gapd_length/2) < (active_length/2)(active_length/2)) gapd_active_count(temp_x,temp_y) = gapd_active_count(temp_x,temp_y) + 1; end end end total_photons = sum(sum(gapd_count));%计算每个细胞中的光子数量（相当于1号要求） active_photons = sum(sum(gapd_active_count));%计算每个细胞中活动区域内的光子数量（相当于2号要求） % 计算PD适用值 result = active_photons/total_photons; final_result = result0.563;%乘以 PD 适用系数获取
基于MATLAB的光纤通信仿真
clc;
clear;
close all;
warning off;
pack;
addpath 'func\'
%是否显示中间结果的波形
%fig1:模拟传感器采集到的数据
%fig2:模拟传感器采集到的数据的净荷值
%fig3:五路传感器数据
%fig4:解封装sync位置
%fig5:显示标准眼图
%fig6:显示噪声干扰眼图
%fig7:显示译码之后眼图
is_show = [0 0 0 0 0 0 0];
Ind = 0;
%参数初始化
EbN0 = [0:0.5:6];
TIMES = 10;
Err_Rate0 = zeros(1,length(EbN0));
Err_Rate0s = zeros(1,TIMES);
Err_Rate1 = zeros(1,length(EbN0));
Err_Rate1s = zeros(1,TIMES);
Err_Rate2 = zeros(1,length(EbN0));
Err_Rate2s = zeros(1,TIMES);
Err_Rate3 = zeros(1,length(EbN0));
Err_Rate3s = zeros(1,TIMES);
Err_Rate4 = zeros(1,length(EbN0));
Err_Rate4s = zeros(1,TIMES);
Err_Rate5 = zeros(1,length(EbN0));
Err_Rate5s = zeros(1,TIMES);
Q1 = zeros(1,length(EbN0));
Q2 = zeros(1,length(EbN0));
Q3 = zeros(1,length(EbN0));
Q4 = zeros(1,length(EbN0));
Q5 = zeros(1,length(EbN0));
Q0 = zeros(1,length(EbN0));
for jj1 = 1:length(EbN0)
jj1
Ind = Ind + 1;
Err_Rate2 = zeros(1,TIMES);
for jj2 = 1:TIMES
jj2
Num_NET = 5;%模拟子网的个数
Lens = 500;%传感器采集数据长度
%传感信息的净荷值相关参数
A = -40;
B = 120;
Delta = 0.2;
Ld = (ceil(log2((B-A)/Delta)));
%传感子网数据帧的封装
Preamble=[0,1,1,1,1,1,1,0];%帧前导码取特殊字符串01111110，表示帧同步，便于传感监控中心判断出帧的起始位置，
%子网ID号，本案为5个子网，所以ID为三bit数据,但是为了具有扩展性，ID用四个bit表示
ID =[0,0,0,0;
0,0,0,1;
0,0,1,0;
0,0,1,1;
0,1,0,0];
%HEC为帧头校验码
HEC = [0,1,0,1];
%包头
Head = [0,0,1,1,0,0,1,1];
%address,地址，随着采集数据，地址逐渐加1
%传感器类型,假设五个子网，每个子网就一种类型
type1 = [0,0,0,1;
0,0,1,0;
0,0,1,1;
0,1,0,0;
0,1,0,1];
%包尾
Trail = [1,1,0,0,1,1,0,0];
Bao_Size= 10;
%最后的复用帧的相关参数
Sync = [0,0,0,0,0,0,1,1,1,0,1,0,1,1,1,1,1,0,0,0,1,0,1,1,0,0,0,0,0,1,0,1,1,1,0,1,0,1,1,1,1,1,0,1,1,1,0,1,0,0,0,0,0,1,0,1,1,1,0,1,0,1,0,0];
Separator = [0,1,0,1];
FCS = [1,0,0,1,1,0,1,1,1,0,0,1,1,0,0,1];
N = 5;
M = 5;
%第一部分，系统发送部分
%第一部分，系统发送部分
%本案共有五个子网，所以模拟五组数据
for nn = 1:Num_NET
%模拟产生传感器采集到的数据,这里，为了模拟不同的码率，将采集的数据长度设置为不同
y = func_sensor_samples(Lens,A,B);
YY{nn} = y;
if is_show(1) == 1
figure(1);
plot(y,'b');
xlabel('times');
ylabel('amplitude')
title('模拟传感器采集到的数据');
pause(0.1);
end
%光纤光栅传感信号的净荷值的二进制转换
Lp = func_value2realvalue(y,A,B,Delta);
if is_show(2) == 1
figure(2);
plot(Lp,'b');
xlabel('times');
ylabel('amplitude')
title('模拟传感器采集到的数据的净荷值');
pause(0.1);
end
%计算转换为二进制
V2 = func_dec2bin(Lp,Ld);
%将二进制转换为串行流
Signal{nn} = (reshape(V2',size(V2,1)*size(V2,2),1))';%通过上面的步骤，我们模拟了实际要发送的二进制码流
end
if is_show(3) == 1
figure(3);
subplot(321);
plot(YY{1});title('第1个子网传感器采集数据波形');
subplot(322);
plot(YY{2});title('第2个子网传感器采集数据波形');
subplot(323);
plot(YY{3});title('第3个子网传感器采集数据波形');
subplot(324);
plot(YY{4});title('第4个子网传感器采集数据波形');
subplot(325);
plot(YY{5});title('第5个子网传感器采集数据波形');
pause(0.1);
end
%Signal即为实际要发送的二进制码流
%产生数据包数据包,这里，我们假设每个数据包中为10个采样数据，即100bit数据
address = 0;
for nn = 1:Num_NET
tmps2 = [];
tmps = [];
address = 0;
for i = 1:length(Signal{nn})/(Bao_Size*Ld)
address = address + 1;
address2 = func_dec2bin_addr(address);
%地址转换为2进制，8bit标识,这样，每个包共128bit
tmps = [Head,address2,type1(nn,:),Signal{nn}((i-1)*Bao_Size*Ld+1:i*Bao_Size*Ld),Trail];
Len_bao = length(tmps);
tmps2 = [tmps2 tmps];
end
Signal2{nn} = tmps2;
end
%由数据帧转换为传感帧
for nn = 1:Num_NET
tmps2 = [];
tmps = [];
for i = 1:length(Signal2{nn})/(N*Len_bao)
tmps = [Preamble,ID(nn,:),HEC,Signal2{nn}((i-1)*Len_bao*N+1:i*Len_bao*N)];
Len_zhen = length(tmps);
tmps2 = [tmps2 tmps];
end
Signal3{nn} = tmps2;
end
%对五个子网的数据进行时分复接
%码率调整
Out1 = func_Rate_sync(Signal3{1});
Out2 = func_Rate_sync(Signal3{2});
Out3 = func_Rate_sync(Signal3{3});
Out4 = func_Rate_sync(Signal3{4});
Out5 = func_Rate_sync(Signal3{5});
%时分复接
Out = func_tdma(Out1,Out2,Out3,Out4,Out5,Len_zhen);
%传感复用帧
tmps2 = [];
tmps = [];
for i = 1:length(Out)/(M*Len_zhen)
tmps = [Sync,Separator,Out((i-1)*M*Len_zhen+1:i*M*Len_zhen),Separator,FCS];
Len_zhen2 = length(tmps);
tmps2 = [tmps2 tmps];
end
Signal4 = tmps2;
%进行编码，编码矩阵为96*192，每次取80,并补充16个0作为虚拟填充符进行编码，共96个数据进行编码
load H;
load G;
Trans_Signal = [];
for i = 1:length(Signal4)/80
Trans_data = [zeros(1,16),Signal4((i-1)*80+1:(i)*80)];
%编码
data_code = mod(Trans_data*G,2);
%BPSK
Trans_BPSK = 2*data_code-1;
Len_code = length(Trans_BPSK);
Trans_Signal = [Trans_Signal,Trans_BPSK];
end
%第二部分，信道，信道部分，由于是光纤传输，且您论文中没有详细说明这个部分内容，所以这里暂时仅考虑高斯白噪声
%第二部分，信道，信道部分，由于是光纤传输，且您论文中没有详细说明这个部分内容，所以这里暂时仅考虑高斯白噪声
sigma = sqrt(1./(10^(EbN0(Ind)/10)*0.5));
Rec_Signal = Trans_Signal + sigma*randn(1,length(Trans_Signal));
%第三部分，监控中心——接收端
%第三部分，监控中心——接收端
%译码
R_Signal = [];
for i = 1:length(Signal4)/80
z_hat = func_Dec(Rec_Signal((i-1)*Len_code+1:i*Len_code),sigma,H,50);
x_hat = z_hat(size(G,2)+1-size(G,1):size(G,2));
R_Signal = [R_Signal,x_hat(17:end)'];%去掉16个填充符
end
if is_show(5) == 1
st = func_eye(Trans_Signal);
eyediagram(st,40)
ylabel('信号幅度');
title('原始信号眼图');
pause(0.1);
end
if is_show(6) == 1
st = func_eye(Rec_Signal);
eyediagram(st,40)
ylabel('信号幅度');
title('接收带噪声干扰信号眼图');
pause(0.1);
end
if is_show(7) == 1
st = func_eye(R_Signal*2-1);
eyediagram(st,40)
ylabel('信号幅度');
title('译码之后信号眼图');
pause(0.1);
end
%解封装
%通过搜索sync来确定每帧的帧头
[R_Signal,Sync_pos] = func_find_sync(R_Signal,Sync);
if is_show(4) == 1
figure(4);
plot(Sync_pos,'b');
hold on;
plot(find(Sync_pos>=40),Sync_pos(find(Sync_pos>=40)),'r*');
hold off;
xlabel('times');
ylabel('amplitude')
title('解封装sync位置');
pause(0.1);
end
%判断当前帧位置能否检测到帧头
if_sync = zeros(1,length(Out)/(M*Len_zhen));%这个变量用来存放是否检测到当前帧，如果为0，则该帧未检测到，直接丢弃，检测下一帧
for i = 1:length(Out)/(M*Len_zhen)
if Sync_pos(1+(i-1)*Len_zhen2) > 40
if_sync(i)=1;
end
end
%解封装
tmps3 = [];
tmps2 = [];
tmps = [];
for i = 1:length(Out)/(M*Len_zhen)
if if_sync(i) == 1
tmps = R_Signal((i-1)*Len_zhen2+1:(i)*Len_zhen2);
else
tmps = zeros(1,Len_zhen2);% 如果该帧没有检测到，那么直接赋值0
end
%解封装
tmps2 = tmps(length(Sync)+length(Separator)+1:M*Len_zhen+length(Sync)+length(Separator));
tmps3 = [tmps3 tmps2];
end
Outs = tmps3;
%数字分接单元
[I1,I2,I3,I4,I5] = func_tdma2(Outs,Len_zhen,M);
%对五个子网的数据进行时分复接
%码率调整
Outr{1} = func_Rate_sync(I1);
Outr{2} = func_Rate_sync(I2);
Outr{3} = func_Rate_sync(I3);
Outr{4} = func_Rate_sync(I4);
Outr{5} = func_Rate_sync(I5);
%将5路传感帧转换为数据帧，进一步拆封
%由数据帧转换为传感帧
for nn = 1:Num_NET
tmps2 = [];
tmps = [];
LL = length(Preamble)+length(ID(nn,:))+length(HEC);
for i = 1:length(Signal2{nn})/(N*Len_bao)
tmps = Outr{nn}((i-1)*Len_zhen+1:i*Len_zhen);
tmps2 = [tmps2 tmps(LL+1:end)];
end
Outr2{nn} = tmps2;
end
%将数据包中数据提取
%Signal即为实际要发送的二进制码流
%产生数据包数据包,这里，我们假设每个数据包中为10个采样数据，即100bit数据
address = 0;
for nn = 1:Num_NET
tmps2 = [];
tmps = [];
LL = length(Head)+length(address2)+length(type1(nn,:));
for i = 1:length(Signal{nn})/(Bao_Size*Ld)
tmps = [Outr2{nn}((i-1)*Len_bao+1:i*Len_bao)];
tmps2 = [tmps2 tmps(LL+1:end-length(Trail))];
end
Outr3{nn} = tmps2;
end
%误码率统计在转换为十进制之前计算
[nberr1,rat1] = biterr(Outr3{1},Signal{1});
[nberr1,rat2] = biterr(Outr3{2},Signal{2});
[nberr1,rat3] = biterr(Outr3{3},Signal{3});
[nberr1,rat4] = biterr(Outr3{4},Signal{4});
[nberr1,rat5] = biterr(Outr3{5},Signal{5});
Err_Rate0(jj2) = (rat1+rat2+rat3+rat4+rat5)/5;
Err_Rate1(jj2) = rat1;
Err_Rate2(jj2) = rat2;
Err_Rate3(jj2) = rat3;
Err_Rate4(jj2) = rat4;
Err_Rate5(jj2) = rat5;
end
Err_Rate0s(Ind) = mean(Err_Rate0);
Err_Rate1s(Ind) = mean(Err_Rate1);
Err_Rate2s(Ind) = mean(Err_Rate2);
Err_Rate3s(Ind) = mean(Err_Rate3);
Err_Rate4s(Ind) = mean(Err_Rate4);
Err_Rate5s(Ind) = mean(Err_Rate5);
Q0(Ind) = sqrt(2)*erfcinv(2*Err_Rate0s(Ind)/10);
Q1(Ind) = sqrt(2)*erfcinv(2*Err_Rate1s(Ind)/10);
Q2(Ind) = sqrt(2)*erfcinv(2*Err_Rate2s(Ind)/10);
Q3(Ind) = sqrt(2)*erfcinv(2*Err_Rate3s(Ind)/10);
Q4(Ind) = sqrt(2)*erfcinv(2*Err_Rate4s(Ind)/10);
Q5(Ind) = sqrt(2)*erfcinv(2*Err_Rate5s(Ind)/10);
disp('over a cycle');
pause(5);
end
figure;
semilogy(EbN0,Err_Rate0s,'b-^');
xlabel('EbN0');
ylabel('Ber');
grid on;
figure;
subplot(151);
semilogy(EbN0,Err_Rate1s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网1误码率');
grid on;
subplot(152);
semilogy(EbN0,Err_Rate2s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网2误码率');
grid on;
subplot(153);
semilogy(EbN0,Err_Rate3s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网3误码率');
grid on;
subplot(154);
semilogy(EbN0,Err_Rate4s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网4误码率');
grid on;
subplot(155);
semilogy(EbN0,Err_Rate5s,'b-^');
xlabel('EbN0');
ylabel('Ber');
title('子网5误码率');
grid on;
figure;
plot(EbN0,Q0,'b-o');
xlabel('EbN0');
ylabel('Q');
title('Q因子');
grid on;