Commit ab17984e authored by alexander.dornheim's avatar alexander.dornheim
Browse files

Started adding CDP_processor. Not yet done.

parent 7986de60
def CDP_processor(config):
# Implementation of the Wirtinger Flow (WF) algorithm presented in the paper
# "Phase Retrieval via Wirtinger Flow: Theory and Algorithms"
# by E. J. Candes, X. Li, and M. Soltanolkotabi
# integrated into the ProxToolbox by
# Russell Luke, September 2016.
# The input data are coded diffraction patterns about a random complex
# valued image.
## Make image
n1 = config['Ny']
n2 = config['Nx'] # for 1D signals, this will be 1
x = randn((n1,n2)) + 1j*randn((n1,n2))
config['truth']=x
config['norm_truth']=norm(x,'fro')
config['truth_dim'] = s.shape
## Make masks and linear sampling operators
L = config['product_space_dimension'] # Number of masks
if n2==1:
Masks = np.random.choice(np.array([1j, -1j, 1, -1]),(n1,L))
elif n1==1:
Masks = np.random.choice(np.array([1j, -1j, 1, -1]),(L,n2))
else:
Masks = zeros((n1,n2,L)); # Storage for L masks, each of dim n1 x n2
# Sample phases: each symbol in alphabet {1, -1, i , -i} has equal prob.
for ll in range(L):
Masks[:,:,ll] = np.random.choice(np.array([1j, -1j, 1, -1]),(n1,n2))
# Sample magnitudes and make masks
temp = rand(Masks.shape);
Masks = Masks * ( (temp <= 0.2)*sqrt(3) + (temp > 0.2)/sqrt(2) )
config['Masks'] = conj(Masks)
# Saving the conjugate of the mask saves on computing the conjugate
# every time the mapping A (below) is applied.
if n2==1:
# Make linear operators; A is forward map and At its scaled adjoint (At(Y)*numel(Y) is the adjoint)
A = @(I) fft(conj(Masks) .* repmat(I,[1 L])); # Input is n x 1 signal, output is n x L array
At = @(Y) mean(Masks .* ifft(Y), 2); # Input is n x L array, output is n x 1 signal
elseif(n1==1)
# Make linear operators; A is forward map and At its scaled adjoint (At(Y)*numel(Y) is the adjoint)
A = @(I) fft(conj(Masks) .* repmat(I,[L 1])); # Input is 1 x n signal, output is L x n array
At = @(Y) mean(Masks .* ifft(Y), 1); # Input is L x n array, output is 1 x n signal
else
A = @(I) fft2(config['Masks .* reshape(repmat(I,[1 L]), size(I,1), size(I,2), L)); # Input is n1 x n2 image, output is n1 x n2 x L array
At = @(Y) mean(Masks .* ifft2(Y), 3); # Input is n1 x n2 X L array, output is n1 x n2 image
end
# Data
Y = abs(A(x));
config['rt_data=Y;
Y=Y.^2;
config['data=Y;
config['norm_data=sum(Y(:))/numel(Y);
normest = sqrt(config['norm_data); # Estimate norm to scale eigenvector
config['norm_rt_data=normest;
## Initialization
npower_iter = config['warmup_iter; # Number of power iterations
z0 = randn(n1,n2); z0 = z0/norm(z0,'fro'); # Initial guess
tic # Power iterations
for tt = 1:npower_iter
z0 = At(Y.*A(z0)); z0 = z0/norm(z0,'fro');
end
Times = toc;
z = normest * z0; # Apply scaling
if(n2==1)
Relerrs = norm(x - exp(-1i*angle(trace(x'*z))) * z, 'fro')/norm(x,'fro');
config['u_0 = repmat(z,1,L);
elseif(n1==1)
Relerrs = norm(x - exp(-1i*angle(trace(z'*x))) * z, 'fro')/norm(x,'fro');
config['u_0 = repmat(z,L,1);
else
Relerrs = norm(x - exp(-1i*angle(trace(x'*z))) * z, 'fro')/norm(x,'fro');
config['u_0=reshape(repmat(z,[1 L]), size(z,1), size(z,2), L);
end
fprintf('Run time of initialization: #.2f\n', Times)
fprintf('Relative error after initialization: #f\n', Relerrs)
fprintf('\n')
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment