壓縮感知的實現（含matlab程式碼）
阿新 • • 發佈：2019-02-15
目錄
---------------------------------------------------------------------------------------------------------------------------------------------
原理參考
（1）劉海英. 基於壓縮感知理論的高光譜影象重建和超分辨成像技術研究[D]. 西安電子科技大學, 2014.
演算法參考
----------------------------------------------------------------------------------------------------------------------------------------------
原理簡介
壓縮感知（Compressive Sensing，CS）。相對於傳統的奈奎斯特取樣定理——要求取樣頻率必須是訊號最高頻率的兩倍或兩倍以上（這就要求訊號是帶限訊號，通常在取樣前使用低通濾波器使訊號帶限），壓縮感知則利用資料的冗餘特性，只採集少量的樣本還原原始資料。
一句話總結我理解的壓縮感知實現方法：以被重建訊號在某個變換域上稀疏作為先驗資訊，用測量矩陣觀測被測訊號，由觀測值結合重建演算法重建出完整的被測訊號。
在具體應用時，我們必須解決 CS 理論的三大關鍵問題：
- 目標訊號的稀疏表示。尋找使得目標訊號 f 變換到其上儘可能稀疏的變換域Ψ ，即訊號稀疏表示問題；
- 測量矩陣的構建。測量矩陣是 CS 理論取樣的實現部分。通過測量矩陣控制的取樣使得目標訊號 f在取樣過程中即被壓縮，同時保證目標訊號所含有效資訊不丟失，能夠由壓縮取樣值還原出目標訊號；
- 重建演算法的設計。重建演算法是從取樣值求解最優化問題尋找到目標訊號最優解。重建演算法的準確性、高效性和穩定性是其設計的關鍵。
對於目標訊號的稀疏表示問題，常見的稀疏基有離散餘弦變換基（DCT）和快速傅立葉變換基（FFT）等。
對於測量矩陣，常見的有高斯隨機矩陣、部分哈達瑪矩陣等。
對於重建演算法，常見的有L1範數、正交匹配追蹤演算法（OMP）等。
對於原理部分，相關文獻、部落格等資源相當多，本文不在這裡贅述，詳情可以參考本文開頭引用內容。
演算法實現
本文分別以稀疏基有離散餘弦變換基（DCT）和快速傅立葉變換基（FFT）做為稀疏基，高斯隨機矩陣、部分哈達瑪矩陣為測量矩陣，L1範數、正交匹配追蹤演算法（OMP）為重建演算法進行壓縮感知演算法實現。
本文以f = cos(2*pi/256*t) + sin(2*pi/128*t)做為原訊號，取原訊號f的20%做為輸入進行壓縮感知重建。
注意：本文main.m中使用了CVX工具箱，CVX工具箱安裝方法參考（CVX工具包（for matlab））
main.m
% 該程式用於驗證壓縮感知理論（包含了L1最小范數求解和OMP求解）
%
%
%
clear all; close all;
%% 產生訊號
choice_transform = 1; % 選擇正交基，1為選擇DCT變換，0為選擇FFT變換
choice_Phi = 0; %選擇測量矩陣，1為部分哈達瑪矩陣，0為高斯隨機矩陣
%-----------------------利用三角函式生成頻域或DCT域離散訊號--------------------------
n = 512;
t = [0: n-1];
f = cos(2*pi/256*t) + sin(2*pi/128*t); % 產生頻域稀疏的訊號
%-------------------------------訊號降取樣率-----------------------
n = length(f);
a = 0.2; % 取原訊號的 a%
m = double(int32(a*n));
%--------------------------------------畫訊號圖--------------------------------------
switch choice_transform
case 1
ft = dct(f);
disp('ft = dct(f)')
case 0
ft = fft(f);
disp('ft = fft(f)')
end
disp(['訊號稀疏度：',num2str(length(find((abs(ft))>0.1)))])
figure('name', 'A Tone Time and Frequency Plot');
subplot(2, 1, 1);
plot(f);
xlabel('Time (s)');
% ylabel('f(t)');
subplot(2, 1, 2);
switch choice_transform
case 1
plot(ft)
disp('plot(ft)')
case 0
plot(abs(ft));
disp('plot(abs(ft))')
end
xlabel('Frequency (Hz)');
% ylabel('DCT(f(t))');
%% 產生感知矩陣和稀疏表示矩陣
%--------------------------利用感知矩陣生成測量值---------------------
switch choice_Phi
case 1
Phi = PartHadamardMtx(m,n); % 感知矩陣（測量矩陣） 部分哈達瑪矩陣
case 0
Phi = sqrt(1/m) * randn(m,n); % 感知矩陣（測量矩陣） 高斯隨機矩陣
end
% Phi = randn(m,n); %randn 生成標準正態分佈的偽隨機數（均值為0，方差為1）
% Phi = rand(m,n); % rand 生成均勻分佈的偽隨機數。分佈在（0~1）之間
f2 = (Phi * f')'; % 通過感知矩陣獲得測量值
% f2 = f(1:2:n);
switch choice_transform
case 1
Psi = dct(eye(n,n)); %離散餘弦變換正交基 程式碼亦可寫為Psi = dctmtx(n);
disp('Psi = dct(eye(n,n));')
case 0
Psi = inv(fft(eye(n,n))); % 傅立葉正變換，頻域稀疏正交基（稀疏表示矩陣）
disp('Psi = inv(fft(eye(n,n)));')
end
A = Phi * Psi; % 恢復矩陣 A = Phi * Psi
%% 重建訊號
%---------------------使用CVX工具求解L1範數最小值-----------------
cvx_begin;
variable x(n) complex;
% variable x(n) ;
minimize( norm(x,1) );
subject to
A*x == f2' ;
cvx_end;
figure;subplot(2,1,2)
switch choice_transform
case 1
plot(real(x));
disp('plot(real(x))')
case 0
plot(abs(x));
disp(' plot(abs(x))')
end
title('Using L1 Norm（Frequency Domain）');
% ylabel('DCT(f(t))'); xlabel('Frequency (Hz)');
switch choice_transform
case 1
sig = dct(real(x));
disp('sig = dct(real(x))')
case 0
sig = real(ifft(full(x)));
disp(' sig = real(ifft(full(x)))')
end
subplot(2,1,1);
plot(f)
hold on;plot(sig);hold off
title('Using L1 Norm (Time Domain)');
% ylabel('f(t)'); xlabel('Time (s)');
legend('Original','Recovery')
%-----------------------------使用OMP演算法重建--------------------------
for K = 1:100
theta = CS_OMP(f2,A,K);
% figure;plot(dct(theta));title(['K=',num2str(K)])
switch choice_transform
case 1
re(K) = norm(f'-(dct(theta)));
case 0
re(K) = norm(f'-real(ifft(full(theta))));
end
end
theta = CS_OMP(f2,A,find(re==min(min(re))));
disp(['最佳稀疏度K=',num2str(find(re==min(min(re))))]);
% theta = CS_OMP(f2,A,10);
figure;subplot(2,1,2);
switch choice_transform
case 1
plot(theta);
disp('plot(theta)')
case 0
plot(abs(theta));
disp('plot(abs(theta))')
end
title(['Using OMP(Frequence Domain) K=',num2str(find(re==min(min(re))))])
switch choice_transform
case 1
sig2 = dct(theta);
disp('sig2 = dct(theta)')
case 0
sig2 = real(ifft(full(theta)));
disp('sig2 = real(ifft(full(theta)))')
end
subplot(2,1,1);plot(f);hold on;
plot(sig2)
hold off;
title(['Using OMP(Time Domain) K=',num2str(find(re==min(min(re))))]);
legend('Original','Recovery')
%%
其中呼叫函式
部分哈達瑪矩陣：PartHadamardMtx.m
function [ Phi ] = PartHadamardMtx( M,N )
%PartHadamardMtx Summary of this function goes here
% Generate part Hadamard matrix
% M -- RowNumber
% N -- ColumnNumber
% Phi -- The part Hadamard matrix
% 來源http://blog.csdn.net/jbb0523/article/details/44700735
%% parameter initialization
%Because the MATLAB function hadamard handles only the cases where n, n/12,
%or n/20 is a power of 2
L_t = max(M,N);%Maybe L_t does not meet requirement of function hadamard
L_t1 = (12 - mod(L_t,12)) + L_t;
L_t2 = (20 - mod(L_t,20)) + L_t;
L_t3 = 2^ceil(log2(L_t));
L = min([L_t1,L_t2,L_t3]);%Get the minimum L
%% Generate part Hadamard matrix
Phi = [];
Phi_t = hadamard(L);
RowIndex = randperm(L);
Phi_t_r = Phi_t(RowIndex(1:M),:);
ColIndex = randperm(L);
Phi = Phi_t_r(:,ColIndex(1:N));
end 正交匹配追蹤演算法：OMP.m
function [ theta ] = CS_OMP( y,A,t )
% 實現壓縮感知OMP演算法
%CS_OMP Summary of this function goes here
%Version: 1.0 written by jbb0523 @2015-04-18
% Detailed explanation goes here
% y = Phi * x
% x = Psi * theta
% y = Phi*Psi * theta
% t 稀疏度
% 令 A = Phi*Psi, 則y=A*theta
% 現在已知y和A，求theta
% 來源：http://blog.csdn.net/jbb0523/article/details/45130793
[y_rows,y_columns] = size(y);
if y_rows<y_columns
y = y';%y should be a column vector
end
[M,N] = size(A);%感測矩陣A為M*N矩陣
theta = zeros(N,1);%用來儲存恢復的theta(列向量)
At = zeros(M,t);%用來迭代過程中儲存A被選擇的列
Pos_theta = zeros(1,t);%用來迭代過程中儲存A被選擇的列序號
r_n = y;%初始化殘差(residual)為y
for ii=1:t%迭代t次，t為輸入引數
product = A'*r_n;%感測矩陣A各列與殘差的內積
[val,pos] = max(abs(product));%找到最大內積絕對值，即與殘差最相關的列
At(:,ii) = A(:,pos);%儲存這一列
Pos_theta(ii) = pos;%儲存這一列的序號
A(:,pos) = zeros(M,1);%清零A的這一列，其實此行可以不要，因為它與殘差正交
%y=At(:,1:ii)*theta，以下求theta的最小二乘解(Least Square)
theta_ls = (At(:,1:ii)'*At(:,1:ii))^(-1)*At(:,1:ii)'*y;%最小二乘解
%At(:,1:ii)*theta_ls是y在At(:,1:ii)列空間上的正交投影
r_n = y - At(:,1:ii)*theta_ls;%更新殘差
end
theta(Pos_theta)=theta_ls;%恢復出的theta
end 
