OMP压缩感知重构仿真

clc;clear%% 1. 时域测试信号生成%产生长度为N=256的稀疏信号，其稀疏度K=23且这23个非零值随机分布于信号256个位置%观测向量y的长度M=80，即采样率M/N=0.3N=256;K=23;M=80;x = zeros(N,1);q = randperm(N);x(q(1:K)) =randn(K,1); %原始信号%% 2. 测量矩阵及观测值获得Phi=randn(M,N); %测量矩阵% 感知矩阵（高斯分布白噪声）M*NmatrixNorm = Phi.'*Phi;matrixNorm = sqrt(diag(matrixNorm)).';Phi = Phi./repmat(matrixNorm, [M,1]); %注意，观测矩阵是要归一化的，因为原子范数要是1！y=Phi*x ; %获得线性测量%% 3.用MP算法重构信号iterations=K; % 算法迭代次数(m>=K)%signal_reconstruct=zeros(1,1); % 近似解矩阵(初始值为空矩阵)r_n=y; % 残差值M*1x_rec=zeros(N,1);for times=1:iterationsfor col=1:N %感知矩阵的所有列向量innerpro(col)=Phi(:,col)'*r_n; %计算余量和感知矩阵每一列的内积end[val,pos]=max(abs(innerpro) ); %找出内积中绝对值最大的元素和它的对应的感知矩阵的列posx_rec(pos)=x_rec(pos)+innerpro(pos); %计算新的近似x_recr_n=r_n-innerpro(pos)*Phi(:,pos); %更新残差endnorm(x_rec-x)/norm(x) % 重构误差subplot(3,1,1);plot(x);title('origin');subplot(3,1,2);plot(x_rec);title('reconstruct');subplot(3,1,3);plot(r_n);title('残差');附录B OMP算法重构一维信号代码clc;clear%% 1. 时域测试信号生成K=7; % 稀疏度(做FFT可以看出来)N=256; % 信号长度M=64; % 测量数(M>=K*log(N/K),至少40,但有出错的概率)f1=50; % 信号频率1f2=100; % 信号频率2f3=200; % 信号频率3f4=400; % 信号频率4fs=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(M,N); % 测量矩阵(高斯分布白噪声)s=Phi*x.'; % 获得线性测量%% 3. 正交匹配追踪法重构信号(本质上是1-范数最优化问题)m=2*K; % 算法迭代次数(m>=K)Psi=fft(eye(N,N))/sqrt(N); % 傅里叶正变换矩阵T=Phi*Psi'; % 恢复矩阵(测量矩阵*正交反变换矩阵)hat_y=zeros(1,N); % 待重构的谱域(变换域)向量Aug_t=[]; % 增量矩阵(初始值为空矩阵)r_n=s; % 残差值for times=1:m; % 迭代次数(有噪声的情况下,该迭代次数为K)for col=1:N; % 恢复矩阵的所有列向量product(col)=abs(T(:,col)'*r_n); % 恢复矩阵的列向量和残差的投影系数(内积值)end[val,pos]=max(product); % 最大投影系数对应的位置Aug_t=[Aug_t,T(:,pos)]; % 矩阵扩充T(:,pos)=zeros(M,1); % 选中的列置零（实质上应该去掉，为了简单我把它置零）aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s; % 最小二乘,使残差最小r_n=s-Aug_t*aug_y; % 残差pos_array(times)=pos; % 纪录最大投影系数的位置endhat_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) % 重构误差附录C OMP算法重构二维信号代码function Wavelet_OMPclc;clear% 读文件X=imread('D:\MATLAB\omp\lena256.bmp');X=double(X);[a,b]=size(X);% 小波变换矩阵生成ww=DWT(a);% 小波变换让图像稀疏化(注意该步骤会耗费时间，但是会增大稀疏度)X1=ww*sparse(X)*ww';X1=full(X1);% 随机矩阵生成M=190;R=randn(M,a);% 测量Y=R*X1;% OMP算法X2=zeros(a,b); % 恢复矩阵for i=1:b % 列循环rec=omp(Y(:,i),R,a);X2(:,i)=rec;end% 原始图像figure(1);imshow(uint8(X));title('原始图像');% 变换图像figure(2);imshow(uint8(X1));title('小波变换后的图像');% 压缩传感恢复的图像figure(3);X3=ww'*sparse(X2)*ww; % 小波反变换X3=full(X3);imshow(uint8(X3));title('恢复的图像');% 误差(PSNR)errorx=sum(sum(abs(X3-X).^2)); % MSE误差psnr=10*log10(255*255/(errorx/a/b)) % PSNR% OMP的函数% s-测量；T-观测矩阵；N-向量大小function hat_y=omp(s,T,N)Size=size(T); % 观测矩阵大小M=Size(1); % 测量hat_y=zeros(1,N); % 待重构的谱域(变换域)向量Aug_t=[]; % 增量矩阵(初始值为空矩阵)r_n=s; % 残差值for times=1:M/4; % 迭代次数(稀疏度是测量的1/4) for col=1:N; % 恢复矩阵的所有列向量product(col)=abs(T(:,col)'*r_n); % 恢复矩阵的列向量和残差的投影系数end[val,pos]=max(product); % 最大投影系数对应的位置Aug_t=[Aug_t,T(:,pos)]; % 矩阵扩充T(:,pos)=zeros(M,1); % 选中的列置零(实质上应该去掉，为了简单我把它置零)aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s; % 最小二乘,使残差最小r_n=s-Aug_t*aug_y; % 残差pos_array(times)=pos; % 纪录最大投影系数的位置if (norm(r_n)<9) % 残差足够小break;endendhat_y(pos_array)=aug_y; % 重构的向量function ww=DWT(N)% 构造正交小波变换矩阵，图像大小N*N，N=2^P，P是整数。

[h,g]= wfilters('sym8','d'); % 分解低通和高通滤波器% N=256; % 矩阵维数(大小为2的整数幂次)L=length(h); % 滤波器长度rank_max=log2(N); % 最大层数rank_min=double(int8(log2(L)))+1; % 最小层数ww=1; % 预处理矩阵% 矩阵构造for jj=rank_min:rank_maxnn=2^jj;% 构造向量p1_0=sparse([h,zeros(1,nn-L)]);p2_0=sparse([g,zeros(1,nn-L)]);% 向量圆周移位for ii=1:nn/2p1(ii,:)=circshift(p1_0',2*(ii-1))';p2(ii,:)=circshift(p2_0',2*(ii-1))';end% 构造正交矩阵w1=[p1;p2];mm=2^rank_max-length(w1);w=sparse([w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)]);ww=ww*w;clear p1;clear p2;end。