%??1-D信號壓縮傳感的實現（正交匹配追蹤法Orthogonal?Matching?Pursuit）
%??測量數M>=K*log（N/K）K是稀疏度N信號長度可以近乎完全重構
%??編程人--香港大學電子工程系?沙威??Email:?wsha@eee.hku.hk
%??編程時間：2008年11月18日
%??文檔下載:?http://www.eee.hku.hk/~wsha/Freecode/freecode.htm?
%??參考文獻：Joel?A.?Tropp?and?Anna?C.?Gilbert?
%??Signal?Recovery?From?Random?Measurements?Via?Orthogonal?Matching
%??Pursuit，IEEE?TRANSACTIONS?ON?INFORMATION?THEORY?VOL.?53?NO.?12
%??DECEMBER?2007.

clc;clear

%%??1.?時域測試信號生成
K=7;??????%??稀疏度（做FFT可以看出來）
N=256;????%??信號長度
M=64;?????%??測量數（M>=K*log（N/K）至少40但有出錯的概率）
f1=50;????%??信號頻率1
f2=100;???%??信號頻率2
f3=200;???%??信號頻率3
f4=400;???%??信號頻率4
fs=800;???%??采樣頻率
ts=1/fs;??%??采樣間隔
Ts=1:N;???%??采樣序列
x=0.3*cos（2*pi*f1*Ts*ts）+0.6*cos（2*pi*f2*Ts*ts）+0.1*cos（2*pi*f3*Ts*ts）+0.9*cos（2*pi*f4*Ts*ts）;??%??完整信號

%%??2.??時域信號壓縮傳感
Phi=randn（MN）;???????????????????????????????????%??測量矩陣（高斯分布白噪聲）
s=Phi*x.‘;????????????????????????????????????????%??獲得線性測量?

%%??3.??正交匹配追蹤法重構信號（本質上是L_1范數最優化問題）
m=2*K;????????????????????????????????????????????%??算法迭代次數（m>=K）
Psi=fft（eye（NN））/sqrt（N）;????????????????????????%??傅里葉正變換矩陣
T=Phi*Psi‘;???????????????????????????????????????%??恢復矩陣（測量矩陣*正交反變換矩陣）

hat_y=zeros（1N）;?????????????????????????????????%??待重構的譜域（變換域）向量?????????????????????
Aug_t=[];?????????????????????????????????????????%??增量矩陣（初始值為空矩陣）
r_n=s;????????????????????????????????????????????%??殘差值

for?times=1:m;????????????????????????????????????%??迭代次數（有噪聲的情況下該迭代次數為K）
????for?col=1:N;??????????????????????????????????%??恢復矩陣的所有列向量
????????product（col）=abs（T（:col）‘*r_n）;??????????%??恢復矩陣的列向量和殘差的投影系數（內積值）?
????end
????[valpos]=max（product）;???????????????????????%??最大投影系數對應的位置
????Aug_t=[Aug_tT（:pos）];???????????????????????%??矩陣擴充
????T（:pos）=zeros（M1）;??????????????????????????%??選中的列置零（實質上應該去掉，為了簡單我把它置零）
????aug_y=（Aug_t‘*Aug_t）^（-1）*Aug_t‘*s;???????????%??最小二乘使殘差最小
????r_n=s-Aug_t*aug_y;????????????????????????????%??殘差
????pos_array（times）=pos;?????????????????????????%??紀錄最大投影系數的位置
end
hat_y（pos_array）=aug_y;???????????????????????????%??重構的譜域向量
hat_x=real（Psi‘*hat_y.‘）;?????????????????????????%??做逆傅里葉變換重構得到時域信號

%%??4.??恢復信號和原始信號對比
figure（1）;
hold?on;
plot（hat_x‘k.-‘）?????????????????????????????????%??重建信號
plot（x‘r‘）???????????????????????????????????????%??原始信號
legend（‘Recovery‘‘Original‘）
norm（hat_x.‘-x）/norm（x）???????????????????????????%??重構誤差