基于EM參數估計的SAGE算法的MATLAB仿真
發布時間：2022/3/15 瀏覽數：5044
這里主要以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');
