30736499CS-based-on-FFT-or-DWT_CS-FFT/DWT_

% 本程序实现图像LENA的压缩传感 % 程序作者：沙威，香港大学电气电子工程学系，wsha@eee.hku.hk % 算法采用正交匹配法，参考文献 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. % 该程序没有经过任何优化 function Wavelet_OMPc clc;clear % 读文件 X=imread('lena256.bmp'); X=double(X); [a,b]=size(X); %# 傅里叶变换矩阵 小波变换矩阵生成 %#ww=DWT(a); %#WW=full(ww); %#正交稀疏基 WW=fft(eye(a,b))/sqrt(a); % 傅里叶变换让图像稀疏化（注意该步骤会耗费时间，但是会增大稀疏度） %X1=ww*sparse(X);%#*ww'; %X1=full(X1);%变换后的含稀疏系数的矩阵 X1=WW'*X; X1=real(X1); %统计稀疏系数的个数 n=0; for i=1:a; for j=1:a; if(X1(i,j)>=1760) n=n+1; end end end K=n % 随机矩阵生成 %M=K*log2(a/K) M=290; R=randn(M,a); % 测量 Y=R*X; %#？为什么测量值用的与含稀疏系数的矩阵的乘积，不是与信号X的乘积？ % OMP算法 X2=zeros(a,b); %#？ 矩阵（%不是恢复矩阵，是含稀疏系数，要重构的矩阵） Z=R*WW; %#传感矩阵 for i=1:b % 列循环 (因为是对矩阵，不是一维信号，所以要进行列循环) rec=omp(Y(:,i),Z,a,K);%#每一列稀疏表示后得到的稀疏向量 X2(:,i)=rec; %将各列稀疏向量组合起来。得到含稀疏系数的矩阵 end % 原始图像 figure(1); imshow(uint8(X)); title('原始图像'); % 变换图像 figure(2); imshow(uint8(X1)); title('小波变换后的图像'); % 压缩传感恢复的图像 figure(3); %X3=ww'*sparse(X2);%*ww; % 小波反变换（将正交稀疏基转置*含稀疏系数的矩阵，得到重构信号） X3=WW*X2; %X3=full(X3); imshow(uint8(X3)); title('恢复的图像'); % 误差(PSNR) errorx=sum(sum(abs(X3-X).^2)); % MSE误差(先计算各对应元素的值的绝对值的平方，然后先将各列的差值加起来，形成一行N列，再将全部差值都加起来。即求所有差值的平方和) psnr=10*log10(255*255/(errorx/a/b)) % PSNR % OMP的函数 % s-测量；T-观测矩阵；N-向量大小 function hat_y=omp(s,T,N,K) %s为测量值向量矩阵的每一列 %？T应为传感矩阵，不是测量矩阵 Size=size(T); % 观测矩阵大小，值为（M，a） M=Size(1); % 测量值数目，为Size的第一个值M hat_y=zeros(1,N); % 待重构的谱域(变换域)向量 Aug_t=[]; % 增量矩阵(初始值为空矩阵) r_n=s; % 残差值 %for times=1:K; 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; end end hat_y(pos_array)=aug_y; % 重构的向量