基于EM参数估计的SAGE算法的MATLAB仿真
发布时间：2022/3/15 浏览数：4271
这里主要以SAGE的基本结构进行算法仿真。整个MATLAB的算法的流程如下所示：
第1步：参数初始化；
第2步：产生发送信号，主要函数是产生论文中的s函数；
第3步：开始SAGE算法循环，转入第4步；
第4步：信号分解；
第5步：延迟估计；
第6步：到达角估计；延迟参数使用前面更新后的值，其余参数使用前一次的迭代值；
第7步：离开角估计；延迟参数和到达角使用前面更新后的值，其余参数使用前一次的迭代值；
第8步：频偏估计；延迟参数、离开角和到达角使用前面更新后的值，其余参数使用前一次的迭代值；
第9步：幅度估计；延迟参数、离开角、频偏和到达角使用前面更新后的值，其余参数使用前一次的迭代值；
第10步：返回第3步，使用更新好的参数再次迭代循环；
其中发送端的基本结构主要根据论文中的图3来实现，图3的基本结构如下所示：
————————————————
部分核心代码如下：
clc;
clear;
close all;
warning off;
addpath 'func\'
%%
%参数初始化
parameter_define;
%变量赋值
Iter = 24; %设置迭代次数
Ker = 100;
Delay = [7*Ker 6*Ker 5*Ker]; %设置延迟，这里，为了提高实际的延迟估计精度，讲每个延迟参数分割k倍
%Theta_out = [2.5133 3.8572 2.3387]/180*pi; %设置离开角
%Theta_in = [4.0666 3.0944 3.9548]/180*pi; %设置进入角
%F_offset = [0.19 0.36 0.38]/180*pi; %设置频偏
%由于这里分辨率较大，论文中也说明了仿真时间需要20多小时，这里
%由于为了验证SAGE算法，所以我们将识别精度设置为可调参数，
%在仿真的时候，适当的降低参数的小数点位数，这样可以加快仿真速度
Theta_out = [0.0439 0.0673 0.0408];
Theta_in = [0.0710 0.0540 0.0610];
F_offset = [0.0023 0.0043 0.0066];
vl = [0.0023 0.0043 0.0066];
Amp = [3.2791 4.5482 3.4650]; %设置幅度
k11 = 0.6;
k12 = 0.9;
k13 = 0.5;
k21 = 1.0;
k22 = 0.7;
k23 = 0.9;
ka = 8.6;
kb = 3.4;
kc = 1.5;
lamda = 0.03; % the length of wave
dt = 0.5*lamda; % the distance between transmit arrays
dr = 0.5*lamda; % the distance between transmit arrays
L = length(Delay); %设置路径数量
Num = 4096; %数据长度
S1 = ones(1,Num);
S2 = ones(1,Num);
S3 = ones(1,Num);
%角度循环范围（理论中为0~2*pi进行扫描，实际中，由于仿真速度和系统内存的限制）
%将范围定为0~0.1（因为原始参数角度均小于0.1，精度为1/1000）
Steps = 0.005;%该值越小，精度越高，误差越小，但仿真时间级数增加
beam = 0:Steps:0.1;%如果电脑够牛逼，那么把0.1修改为2*pi
M = 128; %如果电脑够牛逼，则可改为256,512或者更大
Steps2 = 20000;%频偏精度因子
Delay_set = zeros(L,Iter); %设置延迟迭代变量
Theta_out_set = zeros(L,Iter); %设置离开角迭代变量
Theta_in_set = zeros(L,Iter); %设置进入角迭代变量
F_offset_set = zeros(L,Iter); %设置频偏迭代变量
Amp_set = ones(L,Iter); %设置幅度迭代变量
%%
%产生发送信号
[y,s1,s2,s3,S] = func_S(S1,S2,S3,Amp,F_offset,Delay,Theta_in,Theta_out,Num);
%%
%进行SAGE算法
k = 1;
%初始化
Delay_tmp = zeros(L,1); %设置延迟迭代变量
Theta_out_tmp = zeros(L,1); %设置离开角迭代变量
Theta_in_tmp = zeros(L,1); %设置进入角迭代变量
F_offset_tmp = zeros(L,1); %设置频偏迭代变量
Amp_tmp = ones(L,1); %设置幅度迭代变量
En1 = 0;
En2 = 0;
En3 = 0;
cnt = 1;
%进行迭代
while(k <= Iter - 1)
k
k = k + 1;
%利用SAGE算法开始迭代
%信号分解
%信号分解
%当k = 2的时候，三个路径单独的信号
if k == 2
[s1t,s2t,s3t] = func_Z(S1,S2,S3,F_offset_tmp,Delay_tmp,Theta_in_tmp,Theta_out_tmp,Num);
Z1 = y - s2t - s3t;
Z2 = y - s1t - s3t;
Z3 = y - s1t - s2t;
else
[s1t,s2t,s3t] = func_Z(S1,S2,S3,F_offset_set(:,k-1),Delay_set(:,k-1),Theta_in_set(:,k-1),Theta_out_set(:,k-1),Num);
Z1 = y - s2t - s3t;
Z2 = y - s1t - s3t;
Z3 = y - s1t - s2t;
end
%估算延迟
%估算延迟
%**********************************************************************
if k == 2
Delay_set(:,k) = func_calculate_delay(Z1,Z2,Z3,Num);
else
Delay_set(:,k) = Delay_set(:,k-1);
end
%更新
[s1t,s2t,s3t,a,b,sigma] = func_Z2(S1,S2,S3,F_offset_set(:,k-1),Delay_set(:,k),Theta_in_set(:,k-1),Theta_out_set(:,k-1),Num);
Z1 = y - s2t - s3t;
Z2 = y - s1t - s3t;
Z3 = y - s1t - s2t;
%估算到达角
%估算到达角
%**********************************************************************
count = 0;
for i = [Steps+2*Steps*(k-2):Steps:2*Steps*(k-1)]
count = count + 1;
angle = i;
[tmp1(count),tmp2(count),tmp3(count)] = func_beam_force(Z1,Z2,Z3,Num,angle,M,beam,Steps);
end
Theta_in_set(:,k)= [max(tmp1) max(tmp2) max(tmp3)];
%更新
[s1t,s2t,s3t,a,b,sigma] = func_Z2(S1,S2,S3,F_offset_set(:,k-1),Delay_set(:,k),Theta_in_set(:,k),Theta_out_set(:,k-1),Num);
Z1 = y - s2t - s3t;
Z2 = y - s1t - s3t;
Z3 = y - s1t - s2t;
%估算离开角
%**********************************************************************
count = 0;
for i = [Steps+2*Steps*(k-2):Steps:2*Steps*(k-1)]
count = count + 1;
angle = i;
[tmp1(count),tmp2(count),tmp3(count)] = func_beam_force2(Z1,Z2,Z3,Num,angle,M,beam,Steps);
end
Theta_out_set(:,k)= [max(tmp1) max(tmp2) max(tmp3)];
%更新
[s1t,s2t,s3t,a,b,sigma] = func_Z2(S1,S2,S3,F_offset_set(:,k-1),Delay_set(:,k),Theta_in_set(:,k),Theta_out_set(:,k),Num);
Z1 = y - s2t - s3t;
Z2 = y - s1t - s3t;
Z3 = y - s1t - s2t;
%估算频偏
%**********************************************************************
beamss = [0:1/Steps2:0.01];
[I1tmp,I2tmp,I3tmp,V1tmp,V2tmp,V3tmp,Z1s,Z2s,Z3s] = func_offset_check(s1t,s2t,s3t,a,sigma,b,Steps2,beamss);
F_offset_set(:,k) = [I1tmp/20000 I2tmp/20000 I3tmp/20000];
%更新,做频偏补偿
%Z1s,Z2s,Z3s
%更新
[s1t,s2t,s3t,a,b,sigma] = func_Z2(S1,S2,S3,F_offset_set(:,k),Delay_set(:,k),Theta_in_set(:,k),Theta_out_set(:,k),Num);
Z1 = y - s2t - s3t;
Z2 = y - s1t - s3t;
Z3 = y - s1t - s2t;
%估算幅度
%估算幅度
%**********************************************************************
%幅度值单独由每次迭代得到的四个参数直接计算出,这里做了简化
Amp_set(1,k) = ka*real(s1(end)*exp(-j*2*pi*F_offset_set(1,k)));
Amp_set(2,k) = kb*real(s2(end)*exp(-j*2*pi*F_offset_set(2,k)));
Amp_set(3,k) = kc*real(s3(end)*exp(-j*2*pi*F_offset_set(3,k)));
end
%%
%延迟
figure;
plot(1:Iter,Delay_set(1,:)/Ker,'r-o','LineWidth',2);
hold on;
plot(1:Iter,Delay_set(2,:)/Ker,'b-^','LineWidth',2);
hold on;
plot(1:Iter,Delay_set(3,:)/Ker,'m-*','LineWidth',2);
hold on;
grid on;
axis([1,Iter,0,10]);
legend('路径1','路径2','路径3');
xlabel('迭代次数');
ylabel('延迟收敛值');
%计算最后的误差
disp('(01):延迟——误差百分比%:');
Err11 = 100*abs(Delay_set(1,Iter) - Delay(1))/Delay(1);
Err12 = 100*abs(Delay_set(2,Iter) - Delay(2))/Delay(2);
Err13 = 100*abs(Delay_set(3,Iter) - Delay(3))/Delay(3);
fprintf('%2.4f ',Err11);fprintf('%2.4f ',Err12);fprintf('%2.4f ',Err13);
fprintf('\n---------------------------------\n');
%进入角
figure;
plot(1:Iter,180*Theta_in_set(1,:)/pi,'r-*','LineWidth',2);
hold on;
plot(1:Iter,180*Theta_in_set(2,:)/pi,'b-*','LineWidth',2);
hold on;
plot(1:Iter,180*Theta_in_set(3,:)/pi,'m-*','LineWidth',2);
hold on;
grid on;
legend('路径1','路径2','路径3');
xlabel('迭代次数');
ylabel('进入角收敛值');
%计算最后的误差
disp('(03):进入角——误差百分比%:');
Err31 = 100*abs(Theta_in_set(1,Iter) - Theta_in(1))/Theta_in(1);
Err32 = 100*abs(Theta_in_set(2,Iter) - Theta_in(2))/Theta_in(2);
Err33 = 100*abs(Theta_in_set(3,Iter) - Theta_in(3))/Theta_in(3);
fprintf('%2.4f ',Err31);fprintf('%2.4f ',Err32);fprintf('%2.4f ',Err33);
fprintf('\n---------------------------------\n');
%离开角
figure;
plot(1:Iter,180*Theta_out_set(1,:)/pi,'r-*','LineWidth',2);
hold on;
plot(1:Iter,180*Theta_out_set(2,:)/pi,'b-*','LineWidth',2);
hold on;
plot(1:Iter,180*Theta_out_set(3,:)/pi,'m-*','LineWidth',2);
hold on;
grid on;
legend('路径1','路径2','路径3');
xlabel('迭代次数');
ylabel('离开角收敛值');
%计算最后的误差
disp('(02):离开角——误差百分比%:');
Err21 = 100*abs(Theta_out_set(1,Iter) - Theta_out(1))/Theta_out(1);
Err22 = 100*abs(Theta_out_set(2,Iter) - Theta_out(2))/Theta_out(2);
Err23 = 100*abs(Theta_out_set(3,Iter) - Theta_out(3))/Theta_out(3);
fprintf('%2.4f ',Err21);fprintf('%2.4f ',Err22);fprintf('%2.4f ',Err23);
fprintf('\n---------------------------------\n');
%频偏
figure;
plot(1:Iter,180*F_offset_set(1,:)/pi,'r-*','LineWidth',2);
hold on;
plot(1:Iter,180*F_offset_set(2,:)/pi,'b-*','LineWidth',2);
hold on;
plot(1:Iter,180*F_offset_set(3,:)/pi,'m-*','LineWidth',2);
hold on;
grid on;
legend('路径1','路径2','路径3');
xlabel('迭代次数');
ylabel('频偏收敛值');
%计算最后的误差
disp('(04):频偏——误差百分比%:');
Err41 = 100*abs(F_offset_set(1,Iter) - F_offset(1))/F_offset(1);
Err42 = 100*abs(F_offset_set(2,Iter) - F_offset(2))/F_offset(2);
Err43 = 100*abs(F_offset_set(3,Iter) - F_offset(3))/F_offset(3);
fprintf('%2.4f ',Err41);fprintf('%2.4f ',Err42);fprintf('%2.4f ',Err43);
fprintf('\n---------------------------------\n');
%幅度
figure;
plot(1:Iter,Amp_set(1,:),'r-*','LineWidth',2);
hold on;
plot(1:Iter,Amp_set(2,:),'b-*','LineWidth',2);
hold on;
plot(1:Iter,Amp_set(3,:),'m-*','LineWidth',2);
hold on;
grid on;
legend('路径1','路径2','路径3');
xlabel('迭代次数');
ylabel('幅度收敛值');
%计算最后的误差
disp('(05):幅度——误差百分比%:');
Err51 = 100*abs(Amp_set(1,Iter) - Amp(1))/Amp(1);
Err52 = 100*abs(Amp_set(2,Iter) - Amp(2))/Amp(2);
Err53 = 100*abs(Amp_set(3,Iter) - Amp(3))/Amp(3);
fprintf('%2.4f ',Err51);fprintf('%2.4f ',Err52);fprintf('%2.4f ',Err53);
fprintf('\n---------------------------------\n');
