反演问题 matlab 实现代码

% SCRIPT: iterative_methods % % ------------------------------------------------------------------------- % Comparison of the Iterative Stationary of 1st and 2nd order and CG % methods on a Tikhonov regularization inverse problem: % % observed data: cameraman, blurr unform 9*9, BSNR = 40 dB % % Author: Jose M. Bioucas Dias, 2008 iters = 500; close all; clear alll f = double(imread('lena256.tiff')); %x = double(imread('phantom.tif')); N=length(f); % remove mean mu=mean(f(:)); f=f-mu; % N must be even middle=N/2+1; % blurr matrix B=zeros(N); % Uniform blur lx=4; %blur x-size ly=4; %blurr y-size B((middle-ly):(middle+ly),(middle-lx):(middle+lx))=1; %circularly center B=fftshift(B); %normalize B=B/sum(sum(B)); FB = fft2(B); CFB = conj(FB); FB2 = abs(FB).^2; % define the forward operator A = @(x) real(ifft2(FB.*fft2(x))); AT = @(x) real(ifft2(CFB.*fft2(x))); T = @(x,alpha) real(ifft2(FB2.*fft2(x))) + alpha*x; % convolve g = A(f); % set BSNR BSNR = 40; Pg = var(g(:)); sigma= sqrt((Pg/10^(BSNR/10))); % add noise g=g+ sigma*randn(N); % save figures figure(1); colormap gray; imagesc(f); axis off; figure(2); colormap gray; imagesc(g); axis off; gg = AT(g); % extreme eigenvlues of T tau = 1e-4; x = real(ifft2(CFB./(FB2+tau).*fft2(g))); lam1 = min(FB2(:)) + tau; lamN = max(FB2(:)) + tau; %lam1 = lam1 /10; f_hat = gg; % IM1st clear LI1; beta = 2/(lam1+lamN); for i = 1:iters f_hat = f_hat-beta*(T(f_hat,tau)-gg); LI1(i) = norm(f_hat-x,'fro'); end % IM2nd clear LI2; rho0 = (1-lam1/lamN)/(1+lam1/lamN); alpha = 2/(1+sqrt(1-rho0^2)); beta0 = 2/(lam1+lamN); beta = alpha*beta0; % first iteration f_hat = gg; resid = A(f_hat)-g; f_hat = f_hat-beta0*(T(f_hat,tau)-gg); LI2(1) = norm(f_hat-x,'fro'); f1 = f_hat; f0 = f_hat; for i = 2:iters f_hat = alpha*f1+(1-alpha)*f0-beta*(T(f_hat,tau)-gg); LI2(i) = norm(f_hat-x,'fro'); f0 = f1; f1= f_hat; end % Conjugate gradient method clear LI3; % first iteration f_hat = gg; r = gg -T(f_hat,tau); p = r; for k = 1:iters auxR = r(:)'*r(:); auxT = T(p,tau); a = auxR/(p(:)'*auxT(:)); f_hat = f_hat+a*p; r = r-a*auxT; beta = (r(:)'*r(:))/auxR; p = r + beta*p; LI3(k) = norm(f_hat-x,'fro'); end semilogy(1:length(LI2),LI2, 1:length(LI3),LI3,1:length(LI1),LI1); figure(gcf) title('|f_k-f|^2') legend('2nd order','CG', '1st order')