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

1. Copy each of below codes in different M files.
2. Place all the files in same folder
3. Also note that these codes are not in a particular order. Copy them all and then run the program.
4. 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);

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```