请问这个公式怎么用m​atlab计算并绘图​啊，真心求助。

各位大佬，我在仿真论文里的公式并计算绘图时遇到了，麻烦，论文公式4.3和结果图4.1如图所示：

我自己写的代码如下：

clear;

clc;

syms f;

lamda=1550/10^9;%1550nm

fs=500*10^6;%500MHz

T=400/10^9;%400ns

w=sqrt(log(2))/(log(2)*pi*T*fs);

M=T*fs;

N1=100;

N2=1000;

N3=10000;

CNR=linspace(-40,10,50);%载噪比其实就是信噪比

V_cr1=zeros(1,50);

qianzhixishu1=lamda*fs*w/(2*log(N1*M));%前值系数，就是公式中积分大括号之前的部分

for index=1:50

houzhixishu=(f/w)^2/(1+(CNR(index)/(sqrt(2*pi)*w).*exp(-f^2/(2*w^2)))^-1)^2;%后置系数，其实就是大括号内的积分部分

res=int(houzhixishu,f,-1/2,1/2);

V_cr1(index)=qianzhixishu1*(res^(-1/2));

end

semilogy(CNR,V_cr1);

title('V_c_r_l_b with CNR');

ylabel('V_c_r_l_b(m/s)');

xlabel('CNR(dB)');

结果运行报错：

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

出错 Vcr_CNR_N (第 26 行)

V_cr1(index)=qianzhixishu1*(res^(-1/2));

原因:

错误使用 symengine

Unable to convert expression containing remaining symbolic function calls into double array. Argument must be expression

that evaluates to number.

查了好多网络资料没有结果，真的好心累，谁能帮帮我啊！！！

Answer 1

[qemgafe]

解答如下

clear;

clc;

lamda=1550e-9;%1550nm

fs=500e6;%500MHz

T=400e-9;%400ns

w=sqrt(log(2))/(log(2)*pi*T*fs);

M=T*fs;

N=[100,1000,10000];

for i=1:length(N)

N1=N(i);

CNR1=linspace(-40,10,50);%载噪比其实就是信噪比

CNR=10.^(CNR1/10);

V_cr1=zeros(1,50);

qianzhixishu1=lamda*fs*w/(2*sqrt(N1*M));%前值系数，就是公式中积分大括号之前的部分

for index=1:50

houzhixishu=@(f)(f./w).^2./(1+(CNR(index)./(sqrt(2*pi).*w).*exp(-f.^2./(2.*w.^2)) ).^-1).^2;%后置系数，其实就是大括号内的积分部分

res=integral(houzhixishu,-1/2,1/2);

V_cr1(index)=qianzhixishu1*(res.^(-1/2));

end

semilogy(CNR1,V_cr1);hold on

end

title('V_c_r_l_b with CNR');

ylabel('V_c_r_l_b(m/s)');

xlabel('CNR(dB)');

ylim([1e-3,1e1])

grid on

参数一样的情况下，怀疑文献结果有点问题