% 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')
inverse problem反演问题matlab程序
5星 · 超过95%的资源 需积分: 50 144 浏览量
2010-08-27
17:12:39
上传
评论 15
收藏 56KB RAR 举报
svp_Charles
- 粉丝: 0
- 资源: 2
最新资源
- STM32单片机FPGA毕设电路原理论文报告汽车电动助力转向单片机控制系统设计与试验研究
- STM32单片机FPGA毕设电路原理论文报告气压传感器神经网络算法及单片机实现
- STM32单片机FPGA毕设电路原理论文报告频率的测量在单片机设计中的应用
- 音频转码工具(用于将微信语音 amr 格式转换为 mp3 格式以便在 html5 的 audio 标签中进行播放).zip
- RDK-Web-Performance-Node
- STM32单片机FPGA毕设电路原理论文报告片式电容器浪涌及老化测试系统的设计与实现
- 一个简易的贪吃蛇小游戏(C语言作业).zip
- 国家游戏防沉迷系统接口(php)
- STM32单片机FPGA毕设电路原理论文报告片剂硬度测试仪的液晶显示界面设计
- STM32单片机FPGA毕设电路原理论文报告偏磁式消弧线圈的动态调谐装置
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
- 1
- 2
- 3
- 4
前往页