ABSTRACT
Kalman filter is the recursive data processing algorithm. It generates optimal estimates of desired quantities that are provided as the given set of measurement. Recursive means that property of filter which does not need to store all previous measurement and reprocess all data each time step. In this report we have used Kalman filtering along with the Fractional Fourier and Wavelet transform for the purpose of denoising a real time video. Haar wavelet transform is used . We have used MATLAB as our simulation tool. A real time video is taken as an input. The image array is loaded and then noise is added in the video .After two levels of denoising, our real time video signal gets denoised and a better quality of video signal is achieved. The Kalman filter is a tool that can estimate the variables of a wide range of processes. In the terms of mathematics we can say that a Kalman filter estimates the states of a linear system. It (Kalman filter)does not only works well in practice, but it is theoretically attractive because it can be shown that of all possible filters, it is the one that minimizes the variance of the estimation error. Kalman filters are often implemented in embedded control systems because in order to control a process, we first need an accurate estimate of the process variables.
MATLAB SOURCE CODE
Instructions to run the code
- Copy each of below codes in different M files.
- Place all the files in same folder
- Also note that these codes are not in a particular order. Copy them all and then run the program.
- Run the “FINAL.m” file
Code 1 – Script M File – Final.m
clc clear close all load mri % input the image array % implay(D) orig=D; D=squeeze(D); % preprocessing operations D=im2double(D); D=imnoise(D,'Gaussian',0,0.02); % adding the noise r = size(D,1); % number of rows c = size(D,2); % number of columns f = size(D,3); % number of frames of the video % pause for i=1:f disp(['Frame remaining: ' num2str(f-i)]) xy=D(:,:,i); % ith frame A=wavelet_transform(xy); A1(:,:,i)=A(:,:,1); A2(:,:,i)=A(:,:,2); A3(:,:,i)=A(:,:,3); A4(:,:,i)=A(:,:,4); A5(:,:,i)=A(:,:,5); A6(:,:,i)=A(:,:,6); A7(:,:,i)=A(:,:,7); A8(:,:,i)=A(:,:,8); A9(:,:,i)=A(:,:,9); A10(:,:,i)=A(:,:,10); A11(:,:,i)=A(:,:,11); end [X1,P1,K1]= KALMANFILTER(A1,2); [X2,P2,K2]= KALMANFILTER(A2,2); [X3,P3,K3]= KALMANFILTER(A3,2); [X4,P4,K4]= KALMANFILTER(A4,2); [X5,P5,K5]= KALMANFILTER(A5,2); [X6,P6,K6]= KALMANFILTER(A6,2); [X7,P7,K7]= KALMANFILTER(A7,2); [X8,P8,K8]= KALMANFILTER(A8,2); [X9,P9,K9]= KALMANFILTER(A9,2); [X10,P10,K10]= KALMANFILTER(A10,2); [X11,P11,K11]= KALMANFILTER(A11,2); disp('Press enter to play the original MRI video') pause implay(orig) % displaying the original image disp('Press enter to play the noisy MRI video') pause implay(D) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0') pause implay(X1) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.1') pause implay(X2) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.2') pause implay(X3) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.3') pause implay(X4) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.4') pause implay(X5) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.5') pause implay(X6) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.6') pause implay(X7) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.7') pause implay(X8) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.8') pause implay(X9) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 0.9') pause implay(X10) % displaying the noisy image disp('Press enter to play the de-noised MRI video at fractional angle 1') pause implay(X11) % displaying the noisy image
Code 2 – Function M File -frft.m
function Faf = frft(f, a) % The fast Fractional Fourier Transform % input: f = samples of the signal % a = fractional power % output: Faf = fast Fractional Fourier transform error(nargchk(2, 2, nargin)); f = f(:); N = length(f); shft = rem((0:N-1)+fix(N/2),N)+1; sN = sqrt(N); a = mod(a,4); % do special cases if (a==0), Faf = f; return; end; if (a==2), Faf = flipud(f); return; end; if (a==1), Faf(shft,1) = fft(f(shft))/sN; return; end if (a==3), Faf(shft,1) = ifft(f(shft))*sN; return; end % reduce to interval 0.5 < a < 1.5 if (a>2.0) a = a-2; f = flipud(f); end %%%%%%%%%%%%%%%%%%%%%%%%% if (a>1.5) a = a-1; f(shft,1) = fft(f(shft))/sN; end %%%%%%%%%%%%%%%% if (a<0.5) a = a+1; f(shft,1) = ifft(f(shft))*sN; end %%%%%%%%%%%%%%%%%% % the general case for 0.5 < a < 1.5 alpha = a*pi/2; tana2 = tan(alpha/2); sina = sin(alpha); f = [zeros(N-1,1) ; interp(f) ; zeros(N-1,1)]; % chirp premultiplication chrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2); %both sin and cos terms %chrp = cos(-1*(pi/N*tana2/4*(-2*N+2:2*N-2)'.^2)); %only cos i.e. real terms %chrp = i*sin(-1*(pi/N*tana2/4*(-2*N+2:2*N-2)'.^2)); %only sin i.e. imaginary terms f = chrp.*f; % chirp convolution c = pi/N/sina/4; Faf = fconv(exp(i*c*(-(4*N-4):4*N-4)'.^2),f); %both sin and cos terms %Faf = fconv(cos(c*(-(4*N-4):4*N-4)'.^2),f); %only cos i.e. real terms %Faf = fconv(i*sin(c*(-(4*N-4):4*N-4)'.^2),f); %only sin i.e. imaginary terms Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi); % chirp post multiplication Faf = chrp.*Faf; % normalizing constant Faf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1); %both sin and cos terms %Faf = cos(-1*(1-a)*pi/4)*Faf(N:2:end-N+1); %only cos i.e. real terms %Faf = i*sin(-1*(1-a)*pi/4)*Faf(N:2:end-N+1); %only sin i.e. imaginary terms %%%%%%%%%%%%%%%%%%%%%%%%% function xint=interp(x) % sinc interpolation N = length(x); y = zeros(2*N-1,1); y(1:2:2*N-1) = x; xint = fconv(y(1:2*N-1), sinc([-(2*N-3):(2*N-3)]'/2)); xint = xint(2*N-2:end-2*N+3); %%%%%%%%%%%%%%%%%%%%%%%%% function z = fconv(x,y) % convolution by fft N = length([x(:);y(:)])-1; P = 2^nextpow2(N); z = real(ifft( fft(x,P) .* fft(y,P))); z = z(1:N);
Code 3 – Function M File – haar_decomposition.m
function [a,d]=haar_decomposition(f,v,w) % v=[1/sqrt(2) 1/sqrt(2)]; % w=[1/sqrt(2) -1/sqrt(2)]; m=1:length(f)/2; a=f(2*m-1).*v(1) + f(2*m).*v(2); d=f(2*m-1).*w(1) + f(2*m).*w(2); end
Code 4 – Function M File – haar_reconstruction.m
function recon=haar_reconstruction(a,d,v,w) % v=[1/sqrt(2) 1/sqrt(2)]; % w=[1/sqrt(2) -1/sqrt(2)]; A=[]; D=[]; for i=1:length(a) A=[A a(i)*v]; D=[D d(i)*w]; end recon=A+D; end
Code 5 – Function M File – wavelet_transform.m
function FINAL=wavelet_transform(xy) N=0; FINAL=[]; % decomposition and thresholding % xy=double(xy); for a=0:0.1:1 N=N+1; v1=[1/sqrt(2) 1/sqrt(2)]; % original wavelets w1=[1/sqrt(2) -1/sqrt(2)]; fv=frft(v1,a); %frft of normal scaling function fv=(fv'); fw=([fv(2) -fv(1)]); e=sum(abs(fv).^2); % normalizing the wavelets x=sqrt(1/e); fv1=fv.*x; fw1=([fv1(2) -fv1(1)]); v1=real(fv1); % new wavelets w1=real(fw1); L=[]; H=[]; [r c]=size(xy); for i=1:r f=xy(i,:); [a1,ddd]=haar_decomposition(f,v1,w1); % haar decomposition ddd=thresh_check(ddd); % applying the threshold L(i,:)=a1; H(i,:)=ddd; end [r1,c1]=size(L); LL=[]; LH=[]; HL=[]; HH=[]; for i=1:c1 f1=L(:,i)'; f2=H(:,i)'; [a11,ddd1]=haar_decomposition(f1,v1,w1); [a22,ddd2]=haar_decomposition(f2,v1,w1); ddd1=thresh_check(ddd1); ddd2=thresh_check(ddd2); LL(:,i)=a11'; LH(:,i)=ddd1'; HL(:,i)=a22'; HH(:,i)=ddd2'; end final=[LL LH;HL HH]; % final decomposed image % reconstruction [r2,c2]=size(LL); L=[]; % level 1 reconstruction for j=1:c2 %vertical part1=LL(:,j)'; part2=HL(:,j)'; recon=haar_reconstruction(part1,part2,v1,w1); L=[L recon']; end H=[]; for j=1:c2 %vertical part1=LH(:,j)'; part2=HH(:,j)'; recon=haar_reconstruction(part1,part2,v1,w1); H=[H recon']; end % level 2 reconstruction [r3,c3]=size(L); img=[]; for i=1:r3 % horizontal part1=L(i,:); part2=H(i,:); recon=haar_reconstruction(part1,part2,v1,w1); img(i,:)=recon; end FINAL(:,:,N)=img; end end
Code 6 – Function M File – thresh_check.m
function sig2=thresh_check(sig) M=max(sig); mat=[]; for th=0:0.01:M sig2=sig; for i=1:length(sig) if sig(i)<th && sig(i)>-th sig2(i)=0; end end [PSNR,MSE,MAXERR,L2RAT] = measerr(sig,sig2); mat=[mat MSE]; end [mini,IX]=min(mat); th=0:0.01:M; thresh=th(IX); sig2=sig; for i=1:length(sig) if sig(i)<thresh && sig(i)>-thresh sig2(i)=0; end end end
Code 7 – Function M File – KALMANFILTER.m
function [X,P,K]=KALMANFILTER(video,option) % Usage : % X : the estimated state % P : the estimated error covariance % K : Kalman gain 0<= K <= 1. % The algorithm : % K(n+1) = P_est(n+1) ./ (P_est(n+1)+ Q) . % P_est(n+1) = (I - K).* P_est(n) . % X_est(n+1) = X_est(n) + K .*( Z(n) - H. *X_est(n)) . % option specifies how to estimate the mean state of pixels % % option =1 means that the mean X is calculated on cube [3*3*3] % option =2 means that the mean X is calculated on cube [3*3*2] % Generally the second option preserve edges better % Class of the input Argument "video" : [double] . video(:,:,2:end) = Estimation(video(:,:,2:end),option); [X,P,K] = Perform_Filter(video); %================================Perform_Filter============================ function [X_est,PP,KK] = Perform_Filter(video) % Variables line = size(video,1); % number of rows column = size(video,2); % number of columns time = size(video,3); % number of frames of the video PP = zeros(time,1); % zero column vector of length "time" ie no. of frames KK = zeros(time,1); % zero column vector of length "time" ie no. of frames X_est = double(zeros(size(video))); % empty estimated value which will be converted to the final matrix (video) I = ones(line,column); % matrix of 1s (no of rows)x(no of columns) % Initial conditions . X1 = Initial_Estimation(video(:,:,1)); % estimation of the first frame of "video" X_est(:,:,1) = X1; % putting that value in the actual estimation variable E = X1-video(:,:,1); % difference of the first frame of original video and the estimated video Q = rand(1)/10*ones(line,column); R = rand(1)/10*ones(line,column); P_est = cov(E); % finding the covariance of the difference of the first frame of original video and the estimated video if (line>column) delta = line-column; portion = P_est(1:delta,:); P_est = [P_est;portion]; end if(line<column) P_est = P_est(1:line,:); end K = P_est./(P_est+R); PP(1) = mean(P_est(:)); KK(1) = mean(K(:)); % computation till the first frame % The ALGORITHM for i= 2 : time % working on the rest of the frames of image X_pred = X_est(:,:,i-1); % (i-1)th frame X_pred(isnan(X_pred)) = 0.5; % setting all the nan values to 0.5 X_pred(isinf(X_pred)) = 0.5; % setting all the inf values to 0.5 P_pred = P_est + Q ; % predicted values for ith frame (computed from the estimated value) K = P_pred./(P_pred + R); % Correction = K .*(video(:,:,i)-X_pred); X_est(:,:,i) = X_pred+Correction; X_est(isnan(X_est(:,:,i))) = 0.5; % setting all the nan values to 0.5 X_est(isinf(X_est(:,:,i))) = 0.5; % setting all the nan values to 0.5 P_est = (I-K).*P_pred; PP(i) = mean(P_pred(:)); KK(i) = mean(K(:)); end % eliminating the INF AND NAN in the edges X_est(1,:,:) = X_est(2,:,:); X_est(:,1,:) = X_est(:,2,:); X_est(line,:,:) = X_est(line-1,:,:); X_est(:,column,:)= X_est(:,column-1,:); return %====================================Estimation============================ function Y = Estimation(X,option) % Normalize the data : convert data to Double if ~isa(X,'double') X=double(X)./255; end if max(X(:)) > 1.00 X=X./255.00; end n = size(X,1); % no of rows p = size(X,2); % no of cols m = size(X,3); % no of frames Y = double(zeros(size(X))); % initializing % Estimation % Intial and Final Estimations : X(t_0) an X(t_final) for i=2:n-1 for j=2:p-1 Y(i,j,end) = ( X(i-1,j,m)+X(i+1,j,m)+X(i,j-1,m)+X(i,j+1,m) ) ./ 4; end end if option==1 for i=2:n-1 for j=2:p-1 for k=2:m-1 Y(i,j,k)= (( X(i+1,j-1,k) + X(i+1,j,k)+ X(i+1,j+1,k)+X(i,j-1,k)+... X(i,j+1,k)+X(i-1,j-1,k)+X(i-1,j,k)+... X(i-1,j+1,k) ) + (X(i+1,j-1,k-1) + X(i+1,j,k-1)+ X(i+1,j+1,k-1)+X(i,j-1,k-1)+... X(i,j+1,k-1)+X(i-1,j-1,k-1)+X(i-1,j,k-1)+... X(i-1,j+1,k-1) ) + (X(i+1,j-1,k+1) + X(i+1,j,k+1)+ X(i+1,j+1,k+1)+X(i,j-1,k+1)+... X(i,j+1,k+1)+X(i-1,j-1,k+1)+X(i-1,j,k+1)+... X(i-1,j+1,k+1) )+X(i,j,k) )./27; end end end elseif option==2 for i=2:n-1 for j=2:p-1 for k=2:m Y(i,j,k)=(( X(i+1,j-1,k) + X(i+1,j,k)+ X(i+1,j+1,k)+X(i,j-1,k)+... X(i,j+1,k)+X(i-1,j-1,k)+X(i-1,j,k)+... X(i-1,j+1,k) ) + (X(i+1,j-1,k-1) + X(i+1,j,k-1)+ X(i+1,j+1,k-1)+X(i,j-1,k-1)+... X(i,j+1,k-1)+X(i-1,j-1,k-1)+X(i-1,j,k-1)+... X(i-1,j+1,k-1) +X(i,j,k)) )./18; end end end end % Now all pixels of Y are estimated by 3*3*2 pixels of noisy X % Temprary Stop if (option~=1 &&option~=2) error(' Invalid Option for type of estimation') end return %========================================================================= function Z = Initial_Estimation(X) [n p]=size(X); Z=double(zeros(size(X))); for i=2:n-1 for j=2:p-1 Z(i,j)=X(i-1,j)+X(i+1,j)+X(i,j-1)+X(i,j+1); % 4 connectivity Z(i,j)=Z(i,j)./4; % mean of 4 connectivity pixels end end % Z is now an estimation of X return