1

I want to use CS to reconstruct an image from fewer samples.

I use Gaussian random matrix as measurement matrix. My problem is with Psi matrix which I want to be Haar wavelet coefficients but I don't know how to define it.

I have used DCT and Fourier basis and it worked well. Here is my code with Fourier basis.

Can anyone tell me how to define Psi matrix as haar wavelet transform?

Thanks in advance.

clc
clear all
close all
[fn,fp]=uigetfile({'*.*'});
tic 
A=im2double(rgb2gray(imread([fp,fn])));
figure(1),imshow(A)
xlabel('original')
x=A(:);
n=length(x);
m=1900;
Phi=randn(m,n);   %Measurment Matrix
Psi=fft(eye(n));   %sensing Matrix( or can be dct(eye(n)) )
y=Phi*x;  %compressed signal
Theta=Phi*Psi;
%Initial Guess:  y=Theta*s => s=Theta\y
s2=Theta\y;
%Solution
s1=OMP( Theta, y, 1e-3);
%Reconstruction
x1=Psi*s1;
figure,imshow(reshape(x1,size(A))),xlabel('OMP')
toc
Michael
  • 3,308
  • 5
  • 24
  • 36
Haybert Markarian
  • 47
  • 9

1 Answers1

1

You just need to generate a haar matrix of appropriate dimensions. Consider this MATLAB function:

function [h]=haargen(N)
% Generating Haar Matrix
ih=zeros(N,N); 
h(1,1:N)=ones(1,N)/sqrt(N);
for k=1:N-1 
p=fix(log(k)/log(2)); 
q=k-(2^p); 
k1=2^p; t1=N/k1; 
k2=2^(p+1); t2=N/k2; 
for i=1:t2 
h(k+1,i+q*t1)   = (2^(p/2))/sqrt(N); 
h(k+1,i+q*t1+t2)    =-(2^(p/2))/sqrt(N); 
end

end

Astro
  • 229
  • 2
  • 12
  • Yes I know that but it doesn't work.I mean after reconstruction with OMP algorithm the output is like noisy shape or something like that but when I run the code with DCT or Fourier it works well. – Haybert Markarian Mar 25 '15 at 08:39
  • @HaybertMarkarian This is obvious, because if you see atoms of haar matrix are like square waves i.e. entries are root(2), 1 or -1. And with sparsity constraints your reconstruction will suffer. While for DCT or DFT we have sin/cosine of different frequencies and we are finding a best fourier series kind of representation. – Astro Mar 26 '15 at 08:53
  • @HaybertMarkarian One more thing, with haar matrix you get haar transform coefficients. In wavelets you generally do decomposition in several levels. And the energy packing properties of the Haar transform are not very good. – Astro Mar 26 '15 at 09:03
  • @HaybertMarkarian Alternatively for use of wavelets in CS, we first transforms the signal to obtain coefficient vector, then make it sparse by say thresholding and then sense it with measurement matrix. Finally we try to recover it back. – Astro Mar 26 '15 at 09:11
  • Thanks for your help. Actually I transformed my signal and made it sparse by thresholding and sense it with gaussian matrix. but my problem in this case is with recovering signal. As you know based on Cs we have y=phi*x=phi*psy*s=tetha*s and OMP recovery algorithm which I downloaded from Justin Romberg's site wants me to define psy matrix and I don't know how to define it in wavelet case! cause I only have phi which is gaussian matrix and x which contains sparse wavelet coefficients of my signal. – Haybert Markarian Mar 28 '15 at 09:11
  • See the OMP implementation solves y=As, s being sparse. You can solve two models with this 1) y=Ms, M is measurement matrix 2) y=Mx=MDs, D is dictionary. In 1) A=M and in 2) A=MD. Just call s=OMP(A,y,[],sparsity-K) and recover signal x=Ds. Hope this helps – Astro Mar 28 '15 at 09:25
  • actually theoretically I know Cs based on wavelets but I think I am doing somthing wrong in my matlab code. If it is possible send me your E-mail & I will send you my m-file.Please take a brief look at it, and maybe that way I can find what I am doing wrong with my code. – Haybert Markarian Mar 29 '15 at 12:44
  • Contact me at votrix13@IEEE.org – Astro Mar 30 '15 at 18:18
  • Now you have posted the code. See the problem is full image recovery. Its a general principle: most of real signals follows stationarity in short duration/spatial area only. So try to process you signal in short frames for speech and small patches in case of images. There is no work in cs except one or two which talk about full image recovery. Try to process your image on patch by patch bases e.g., 8x8. You will get desired results. More over with omp you can do batch peocessing for all patches simulatenously. Hope this helps – Astro Mar 30 '15 at 18:25