clc;clear;
close?all;
load??MIG25.MAT
X_fft_rng=fft（X）;
X_fft_rng=fftshift（X_fft_rng1）;
data=X_fft_rng（:385:512）;

K=13;?????%??稀疏度
N=128;????%??信號長度
M=64;????%??測量數（M>=K*log（N/K）但有出錯的概率）

%%?2-D?FFT經典成像方法
%?對完成相位補償的信號進行方位維fft成像
comed_echo=data;%沒有采樣損失的原始數據
temp=fft（comed_echo[]2）;%方位維FFT
img_cls=fftshift（temp2）;%方位維翻折

%?將完整數據的成像結果繪制成圖
figure（1）;
imagesc（abs（img_cls））
%?imshow（1-flipud（log10（20*（imm/max（max（imm））））））;
axis?on;
xlabel（‘方位向‘）;
ylabel（‘距離向‘）;
title（‘完整數據的2D成像結果‘）

%?Phi?=?randn（MN）;
%?Psi?=?fft（eye（NN））/sqrt（N）;
%?T?=?Phi*Psi‘;
data=awgn（data10）;
img_CS?=?zeros（64N）;
Error?=?zeros（164）;%64個方位維的重構誤差
for?iCol=1:64
????x=data（iCol:）;
????x?=?x.‘;
????%?時域信號壓縮傳感
????Phi=randn（MN）;???????????????????????????????????%??測量矩陣（高斯分布白噪聲）
????y=Phi*x;??????????????????????????????????????????%??獲得線性測量
????%?正交匹配追蹤法重構信號（本質上是L_1范數最優化問題）
????m=K;??????????????????????????????????????????????%??算法迭代次數（m>=K）
????Psi=fft（eye（NN））/sqrt（N）;????????????????????????%??傅里葉正變換矩陣
????T=Phi*Psi‘;???????????????????????????????????????%??恢復矩陣（測量矩陣*正交反變換矩陣）
????s_hat=zeros（N1）;?????????????????????????????????%??待重構的譜域（變換域）向量
????Aug_t=[];?????????????????????????????????????????%??增量矩陣（初始值為空矩陣）
????res_n=y;??????????????????????????????????????????%??殘差值
????for?iter=1:m??????????????????????????????????????%??迭代次數（有噪聲的情況下該迭代次數為K）
????????product=abs（T‘*res_n）;????????????????????????%??恢復矩陣的列向量和殘差的投影系數（內積值）
????????[valpos]=max（product）;???????????????????????%??最大投影系數對應的位置
????????Aug_t=[Aug_tT（:pos）];???????????????????????%??矩陣擴充
????????T（:pos）=zeros（M1）;??????????????????????????%??選中的列置零（實質上應該去掉，為了簡單我把它置零）
????????aug_y=（Aug_t‘*Aug_t）^（-1）*Aug_t‘*y;???????????%??最小二乘使殘差最小
????????res_n=y-Aug_t*aug_y;??????????????????????????%??殘差
????????pos_array（iter）=pos;??????????????????????????%??紀錄最大投影系數的位置
????end
????s_hat（pos_array）=aug_y;???????????????????????????%??重構的譜域向量
????img_CS（iCol:）=sqrt（N）*s_hat;
????Error（iCol）?=?norm（s_hat-img_cls（iCol:））/norm（img_cls（iCol:））;
end
img_CS=fftshift（img_CS2）;
figure（2）;
imagesc（abs（img_CS））
title（‘50%稀疏成像結果‘）;
axis?on;
xlabel（‘方位向‘）;
ylabel（‘距離向‘）;

%?誤差分析
norm（img_CS）/norm（img_cls）
figure（3）

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件??????524310??2020-04-16?15:22??MIG25.MAT
?????文件????????2578??2020-04-17?16:43??MIG25CS_OMP1.m