clc;
clearvars;
close all;
% 读文件
X=imread('mandrill256.bmp');
tic;
X=double(X);
[m,n]=size(X);
% % 小波变换矩阵生成
[LL1, LH1, HL1, HH1] = dwt2(X, 'haar');
[LL2, LH2, HL2, HH2] = dwt2(LL1, 'haar');
% [LL3, LH3, HL3, HH3] = dwt2(LL2, 'haar');
% [LL4, LH4, HL4, HH4] = dwt2(LL3, 'haar');
% [LL5, LH5, HL5, HH5] = dwt2(LL4, 'haar');
% LL4 = [LL5, LH5; HL5, HH5];
% LL3 = [LL4, LH4; HL4, HH4];
% LL2=[LL3, LH3; HL3, HH3];
LL1=[LL2, LH2; HL2, HH2];
X1=[LL1, LH1; HL1, HH1];
% 随机矩阵生成
M=m/4;
H = hadamard(m);
R = H(1:M,:);
% 将X1更加稀疏化处理
X1_show=X1;
key_value=find_value_at_a_percent(X1,(M/m)/4);
X1(abs(X1)<key_value)=0;
% 测量值
Y=R*X1;
% OMP算法
% 恢复矩阵
X2=zeros(m,n);
% 按列循环
for i=1:n
% 通过OMP，返回每一列信号对应的恢复值（小波域）
% rec=omp(Y(:,i),R,m);
rec=omp(R,Y(:,i),m);
% rec=ompOptimized(Y(:,i),R,m);
% rec=orthogonal_matching_pursuit(Y(:,i),R,m);
% 恢复值矩阵，用于反变换
X2(:,i)=rec;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 执行逆小波变换
[LL1, LH1, HL1, HH1] = partitionMatrix(X2);
[LL2, LH2, HL2, HH2] = partitionMatrix(LL1);
% [LL3, LH3, HL3, HH3] = partitionMatrix(LL2);
% [LL4, LH4, HL4, HH4] = partitionMatrix(LL3);
% [LL5, LH5, HL5, HH5] = partitionMatrix(LL4);
% 第5级逆变换
% LL4=idwt2(LL5, LH5, HL5, HH5, 'haar');
% 第四级逆变换
% LL3 = idwt2(LL4, LH4, HL4, HH4, 'haar');
% 第三级逆变换
% LL2= idwt2(LL3, LH3, HL3, HH3, 'haar');
% 第二级逆变换
LL1 = idwt2(LL2, LH2, HL2, HH2, 'haar');
X3 = idwt2(LL1, LH1, HL1, HH1, 'haar');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 显示部分
subplot(2,2,1),imshow(uint8(X));
title('原始图像');
subplot(2,2,2),imshow(X1,[]);
title("X1");
subplot(2,2,3),imshow(X2,[]);
title("X2");
subplot(2,2,4),imshow(uint8(X3),[]);
title('恢复的图像');
% 误差(PSNR)
mypsnr=psnr(uint8(X3),uint8(X));
disp(['psnr:' num2str(mypsnr)]);
toc;
% OMP的函数
% s-测量；T-观测矩阵；N-向量大小
function [X1, X2, X3, X4] = partitionMatrix(X)
% 获取矩阵X的大小
[m, n] = size(X);
% 将矩阵X划分为四等份
X1 = X(1:m/2, 1:n/2); % 左上角子矩阵
X2 = X(1:m/2, n/2+1:end); % 右上角子矩阵
X3 = X(m/2+1:end, 1:n/2); % 左下角子矩阵
X4 = X(m/2+1:end, n/2+1:end); % 右下角子矩阵
end
function value_at_a_percent = find_value_at_a_percent(matrix,aPercent)
% 对矩阵进行降序排序
sorted_matrix = sort(matrix(:), 'descend');
% 计算15%位置的索引
num_elements = numel(sorted_matrix);
index_a_percent = ceil(aPercent * num_elements);
% 获取15%位置的值
value_at_a_percent = sorted_matrix(index_a_percent);
end
% OMP算法实现
function x_hat = omp(A, y, M)
N = size(A, 2);
x_hat = zeros(N, 1);
residual = y;
selected_indices = [];
tolerance = 1e-6; % 设置残差阈值
for i = 1:M % 为避免过拟合，迭代次数不能超过观测数量
% 计算测量残差的投影
projections = abs(A' * residual);
% 选择最大投影对应的索引
[~, index] = max(projections);
% 添加选定的原子索引
selected_indices = [selected_indices, index];
% 更新估计的稀疏信号
x_hat(selected_indices) = A(:, selected_indices) \ y;
% 更新残差
residual = y - A * x_hat;
% 检查残差是否足够小//停止条件
if norm(residual) < tolerance
break;
end
end
end
运行结果
dct实现
clc;
clearvars;
close all;
% 读文件
X=imread('mandrill256.bmp');
tic;
X=double(X);
[m,n]=size(X);
% % 小波变换矩阵生成
X1=dct2(X);
% 随机矩阵生成
M=m/4;
H = hadamard(m);
R = H(1:M,:);
% 将X1更加稀疏化处理
X1_show=X1;
key_value=find_value_at_a_percent(X1,(M/m)/4);
X1(abs(X1)<key_value)=0;
% 测量值
Y=R*X1;
% OMP算法
% 恢复矩阵
X2=zeros(m,n);
% 按列循环
for i=1:n
% 通过OMP，返回每一列信号对应的恢复值（小波域）
% rec=omp(Y(:,i),R,m);
rec=omp(R,Y(:,i),m);
% rec=ompOptimized(Y(:,i),R,m);
% rec=orthogonal_matching_pursuit(Y(:,i),R,m);
% 恢复值矩阵，用于反变换
X2(:,i)=rec;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 逆向变换
X3 = idct2(X2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 显示部分
subplot(2,2,1),imshow(uint8(X));
title('原始图像');
subplot(2,2,2),imshow(X1,[]);
title("X1");
subplot(2,2,3),imshow(X2,[]);
title("X2");
subplot(2,2,4),imshow(uint8(X3),[]);
title('恢复的图像');
% 误差(PSNR)
mypsnr=psnr(uint8(X3),uint8(X));
disp(['psnr:' num2str(mypsnr)]);
toc;
% OMP的函数
% s-测量；T-观测矩阵；N-向量大小
function value_at_a_percent = find_value_at_a_percent(matrix,aPercent)
% 对矩阵进行降序排序
sorted_matrix = sort(matrix(:), 'descend');
% 计算15%位置的索引
num_elements = numel(sorted_matrix);
index_a_percent = ceil(aPercent * num_elements);
% 获取15%位置的值
value_at_a_percent = sorted_matrix(index_a_percent);
end
% OMP算法实现
function x_hat = omp(A, y, M)
N = size(A, 2);
x_hat = zeros(N, 1);
residual = y;
selected_indices = [];
tolerance = 1e-6; % 设置残差阈值
for i = 1:M % 为避免过拟合，迭代次数不能超过观测数量
% 计算测量残差的投影
projections = abs(A' * residual);
% 选择最大投影对应的索引
[~, index] = max(projections);
% 添加选定的原子索引
selected_indices = [selected_indices, index];
% 更新估计的稀疏信号
x_hat(selected_indices) = A(:, selected_indices) \ y;
% 更新残差
residual = y - A * x_hat;
% 检查残差是否足够小//停止条件
if norm(residual) < tolerance
break;
end
end
end
输出
参考文献：
https://blog.csdn.net/qq_43095255/article/details/135087322
