Noise removal in a signal using newly designed wavelets and kalman filtering

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);
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

 

Write Your Comments

Your email address will not be published. Required fields are marked *

Recent Posts

Tags

ad-hoc networks AODV boundary detection process classification clustering clustering algorithm Colour Information computer vision Decryption Encryption EZRP ICM (Iterated Conditional Modes) image denoising image enhancement IMAGE PROCESSING image segmentation Imaging and image processing MANET Markov Random Fields neutrosophic logic optical network proposed method PSNR QLab system region growing Robert’s operator Seed point selection segmentation semi-automatic algorithm Shadow Detection shadow removal wall motion wireless communication Wireless network wireless networks Wireless Sensor Network wireless sensor networks ZRP