tuxianghuifu.rar_共轭梯度_共轭梯度 图像_共轭梯度法_图像 谱

% Set parameters. im_name = 'Pirate.tif'; % Name of clean image. noise_std = 3.0; % Standard deviation of image noise. gauss_std = 4.0; % St. dev. of Gaussian PSF. tau =0.432; % Fudge factor in residual bound. rho=0.325; tic; % Load the clean image. Xc = double(imread('LENA.bmp')); %Xc=double(Xc(128:256 ,128:256)); [p,q] = size(Xc); n=p*q; % Construct a Gaussian point spread function. wx = exp(-((1:p)-ceil(p/2)).^2/(2*gauss_std^2))'; wy = exp(-((1:q)-ceil(q/2)).^2/(2*gauss_std^2))'; PSF = wx*wy'; PSF = PSF/sum(sum(PSF)); % Normalized PSF. % Compute the DCT spectrum of the doubly symmetric PSF. e1 = zeros(p,q); e1(1,1) = 1.0; A = dctshift(PSF,[ceil(p/2),ceil(q/2)]); S = dcts2(A)./dcts2(e1); % Blur the clean image and add Gaussian noise; make sure the pixels of the % noisy image are in the range 0,...,255. B = idcts2(dcts2(Xc).*S) + randn(p,q)*noise_std; B(B<0) = 0; B(B>255) = 255; bhat=dcts2(B); bhat=bhat(:); s=S(:); X_iter=bhat; alpha=2; sigma=0.8; k=1; conv=0; TimeCost(1)=0; t0 = cputime; for k=1:4 %while 1 X_old=X_iter; D=s.*(X_old)-bhat; gk=(conj(s).*D+2*tau*X_old); if(norm(gk)<0.1) end if(k==1) p=-gk; else s=(1-sigma)*norm(gk)^2; t=p0'*yk; u=(1-sigma)*gk'*p0; v=p0'*yk; beta=s/t; theta=1+ u/v; p= -theta*gk+ beta*p0; end for i=1:10 df1=gk'*p; fun_val=0.5*norm(s.*(X_old)-bhat)^2+tau*norm((X_old))^2; fun=0.5*norm(s.*(X_old+alpha*p)-bhat)^2+tau*norm(X_old+alpha*p)^2; if fun_val-fun >= -alpha*rho*df1%满足条件1 Y=X_old+alpha*p; D=s.*Y-bhat; g_new=conj(s).*D+2*tau*Y; dff=p'*g_new; if dff >= sigma*df1 %满足条件2 conv=1; break; else alpha = 0.5*alpha; end else alpha = 0.5*alpha; end end if conv==0 disp('未找到满足WOLFE条件的步长，将返回默认步长 alpha=1'); end X_next=X_old+alpha*p; sk=X_next-X_old; D=s.*(X_next)-bhat; g_new=conj(s).*D+2*tau.*X_next; yk=g_new-gk; X_iter=X_next; p0=p; end TimeCost(k+1)=cputime-t0; X_hat=reshape(X_next,size(B)); X=idcts2(X_hat); rmse=norm(Xc-X)/norm(Xc) toc; figure(1), clf, colormap(gray) imagesc(Xc); axis image off title('Original image') figure(2), clf, colormap(gray); imagesc(B); axis image off title('Noisy and blurred image') figure(3), clf, colormap(gray); imagesc(X); axis image off title('deblurred image')