Upsampling in Frequency Domain

1.      Upsampling in Spatial Domain

1.1  Upsampling using MATLAB built-in function

1.2  Upsampling in 1D

1.3  Upsampling in 2D or image matrix

2.      Upsampling in Frequency Domain

2.1  Upsampling a 1D signal

2.2  Upsampling a image matrix

In Frequency domain, upsampling means nothing but the padding of zeros at the end of high frequency components on both sides of the signal.

STEPS TO PERFORM:

1.      Read an image

2.      Obtain the ratio to upsample

3.      Perform Fast Fourier Transform

4.      Shift the Low frequency components to the centre and High frequency components outside.

5.      Add zeros on both the sides of the image

6.      Shift the High frequency components to the centre and Low frequency components to the exterior (Inverse of fftshift)

7.      Perform Inverse Fast Fourier Transform

8.      Display the Upsampled Image

MATLAB CODE:

%UPSAMPLING IN FREQUENCY DOMAIN

%1D UPSAMPLING

FS = 100;

t  = 0:(1/FS):1;

A = 10*sin(2*pi*5*t);

figure,plot(t,A);

 

%FOURIER DOMAIN

FT = fft(A);

fq =linspace(-1/FS,1/FS,101);

figure,plot(fq,abs(FT));title('Before FFTSHIFT');

FT_s = fftshift(FT);

figure,plot(fq,abs(FT_s));title('After FFTSHIFT');

pad_zero = padarray(FT_s,[0 50]);

fq =linspace(-1/FS,1/FS,201);

figure,plot(fq,abs(pad_zero));title('After PADDING WITH ZEROS');

%INVERSE FOURIER TRANSFORM

IFT = ifft(ifftshift(pad_zero));

%UPSAMPLED SIGNAL

t1 = linspace(0,1,201);

figure,plot(t1,(IFT*2),'r',t,A,'g');

EXPLANATION:

Amplitude of the input original signal is 10 and the frequency is 5.

Similarly, the amplitude of the upsampled signal is 10 and the frequency is 5. The number of samples used to plot the signal is increased in the later case.

IMAGE UPSAMPLING IN FOURIER DOMAIN

MATLAB CODE:

%READ AN INPUT IMAGE

A = imread('cameraman.tif');

%RATIO

RatioM = 3;

RatioN = 2;

%UPSAMPLING OVER EACH ROW

mnrow = round(size(A,2)*(RatioM-1)/2);

% 1D FFT ON EACH ROW

row_fft = fft(A,[],2);

%PAD WITH ZEROS ON BOTH SIDES OF EACH ROW

pad_row = padarray(fftshift(row_fft,2),[0 mnrow]);

 

Logarthmic scale

%UPSAMPLING OVER EACH COLUMN

mncol = round(size(A,1)*(RatioN-1)/2);

% 1D FFT ON EACH COLUMN

col_fft = fft(pad_row,[],1);

%PAD WITH ZEROS ON BOTH SIDES OF EACH COLUMN

pad_col = padarray(fftshift(col_fft,1),[mncol 0]);

 

Logarthmic scale

%PERFORM 1D IFFT ON EACH COLUMN

ifft_col = ifft(ifftshift(pad_col,1),[],1);

%PERFORM 1D IFFT ON EACH ROW

ifft_col_row = ifft(ifftshift(ifft_col,2),[],2);

%DISPLAY THE IMAGE AFTER UPSAMPLING

res = abs(ifft_col_row);

res = uint8(res*(numel(res)/numel(A)));

figure,imagesc(res);

SIMPLE VERSION :

A = imread('cameraman.tif');

%RATIO

RatioM = 3;

RatioN = 2;

mnrow = round(size(A,2)*(RatioM-1)/2);

mncol = round(size(A,1)*(RatioN-1)/2);

%FFT ON 2D MATRIX

FT = fftshift(fft2(A));

%PADDING WITH ZEROS

pad_rc = padarray(FT,[mncol mnrow]);

%INVERSE FOURIER TRANSFORM

IFT = ifft2(ifftshift(pad_rc));

Img = uint8(abs(IFT)*(numel(IFT)/numel(A)));

figure,imagesc(Img);

Also Check : UPSAMPLING IN SPATIAL DOMAIN:

Convolution in MATLAB

                 Let us try to understand convolution by performing spatial averaging on a matrix without using MATLAB built in function ‘conv2()’.

%CONVOLUTION IN MATLAB

%INPUT  MATRIX

A = zeros(5);

A(:) = 1:25;

%KERNEL

avg3 = ones(3)/9;

%CONVOLUTION

Result = conv2(A,avg3,'same');

display(Result);

Steps to be performed:

NOTE :

To define a kernel for spatial averaging, fill the kernel with ones and divide it by the number of elements in it.

For instance, consider kernel of size 4x4 , fill the  matrix with ones and divide it by 16. i.e the total number of elements in the matrix.

MATLAB CODE:

%INPUT  MATRIX

A = zeros(5);

A(:) = 1:25;

%KERNEL

avg3 = ones(3)/9;

%PAD THE MATRIX WITH ZEROS

B = padarray(A,[1 1]);

% PRE-ALLOCATE THE MATRIX

Output = zeros([size(A,1) size(A,2)]);

%PERFORM COONVOLUTION

for i = 1:size(B,1)-2

    for j = 1:size(B,2)-2

        Temp = B(i:i+2,j:j+2).*avg3;

        Output(i,j) = sum(Temp(:));

    end

 end

display(Output);

Color Histogram Equalization - MATLAB CODE

                          Histogram Equalization can be considered as redistribution of the intensity of the image. Color histogram equalization can be achieved by converting a color image into HSV/HSI image and enhancing the Intensity while preserving hue and saturation components. 

                          However, performing histogram equalization on components of R,G and B independently will not  enhance the image. At the end of this post, check the histogram of before and after histogram equalization of an image which is obtained by performing histogram equalization on the components(R,G and B) independently.

Steps to be performed:

1.      Convert RGB image into HSI Image.
https://www.imageeprocessing.com/2013/05/converting-rgb-image-to-hsi.html

2.      Obtain the ‘Intensity Matrix’ from the HSI Image matrix

3.      Perform Histogram Equalization on the intensity Matrix
https://www.imageeprocessing.com/2011/04/matlab-code-histogram-equalization.html

4.      Update the Intensity Matrix from the HSI Image matrix with the histogram equalized Intensity matrix

5.      Convert HSI Image back to RGB Image
https://www.imageeprocessing.com/2013/06/convert-hsi-image-to-rgb-image.html

        

MATLAB CODE:

%COLOR HISTOGRAM EQUALIZATION

%READ THE INPUT IMAGE

I = imread('football.jpg');

%CONVERT THE RGB IMAGE INTO HSV IMAGE FORMAT

HSV = rgb2hsv(I);

%https://www.imageeprocessing.com/2013/05/converting-rgb-image-to-hsi.html

%PERFORM HISTOGRAM EQUALIZATION ON INTENSITY COMPONENT

Heq = histeq(HSV(:,:,3));

%https://www.imageeprocessing.com/2011/04/matlab-code-histogram-equalization.html

HSV_mod = HSV;

HSV_mod(:,:,3) = Heq;

RGB = hsv2rgb(HSV_mod);

%https://www.imageeprocessing.com/2013/06/convert-hsi-image-to-rgb-image.html

figure,subplot(1,2,1),imshow(I);title('Before Histogram Equalization');

       subplot(1,2,2),imshow(RGB);title('After Histogram Equalization');

EXPLANATION:

RGB image matrix is converted into HSI(Hue ,Saturation and Intensity) format and histogram equalization is applied only on the Intensity matrix . The Hue and Saturation matrix remains the same. The updated HSI image matrix is converted back to RGB image matrix.

%DISPLAY THE HISTOGRAM OF THE ORIGINAL AND THE EQUALIZED IMAGE

HIST_IN = zeros([256 3]);

HIST_OUT = zeros([256 3]);

%http://angeljohnsy.blogspot.com/2011/06/histogram-of-image.html

%HISTOGRAM OF THE RED,GREEN AND BLUE COMPONENTS

HIST_IN(:,1) = imhist(I(:,:,1),256); %RED

HIST_IN(:,2) = imhist(I(:,:,2),256); %GREEN

HIST_IN(:,3) = imhist(I(:,:,3),256); %BLUE

HIST_OUT(:,1) = imhist(RGB(:,:,1),256); %RED

HIST_OUT(:,2) = imhist(RGB(:,:,2),256); %GREEN

HIST_OUT(:,3) = imhist(RGB(:,:,3),256); %BLUE

mymap=[1 0 0; 0.2 1 0; 0 0.2 1];

figure,subplot(1,2,1),bar(HIST_IN);colormap(mymap);legend('RED CHANNEL','GREEN CHANNEL','BLUE CHANNEL');title('Before Applying Histogram Equalization');

       subplot(1,2,2),bar(HIST_OUT);colormap(mymap);legend('RED CHANNEL','GREEN CHANNEL','BLUE CHANNEL');title('After Applying Histogram Equalization');

EXPLANATION:

Obtain the histogram of each component (Red,Green and Blue) independently.

Define the colormap ‘mymap’ with three colors namely Red, Green and Blue.

Display the histograms of the components before and after histogram equalization.

 NOTE:

Histogram of the above image by processing the components independently gives bad result.

 

Histogram Equalization applied on individual components

Fast Fourier Transform on 2 Dimensional matrix using MATLAB

Fast Fourier transformation on a 2D matrix can be performed using the MATLAB built in function ‘fft2()’.

Fourier transform is one of the various mathematical transformations known which is used to transform signals from time domain to frequency domain.
The main advantage of this transformation is it makes life easier for many problems when we deal a signal in frequency domain rather than time domain.

Example:

%FOURIER TRANSFORM ON A MATRIX

A  = zeros(5);

A ( : ) = 1:25;

display(A);

F_FFT  = fft2(A);

display(F_FFT);

%INVERSE FOURIER TRANSFORM

I_FFT = ifft2(F_FFT);

display(abs(I_FFT));

NOTE on ABSOLUTE VALUE:

When we use FFT2() or FFT(),the result we obtain in the frequency domain is of complex data type.

i.e It contains both the real  as well as the imaginary part.

Let   A=10+5i

A is a complex number as it contains both real and imaginary part.In this particular case ‘10’ is the real part and ‘5’ is the imaginary part.

abs(A) = 11.1803 is the absolute (also called modulus in few books or notations) value of A which is nothing but the magnitude. It can be arrived by using the below mentioned formula:

abs(A) = sqrt(real part^2+imaginary part^2).

              = sqrt(10^2+5^2)

              = sqrt(125)

             = 11.1803 (approx)

Let’s try to understand how the Fourier transform on 2 dimensional data works with a simple example.

This method will be helpful to understand the up sampling and down sampling in both spatial and frequency domain.

1.       Consider a matrix A

2.       Perform 1 D Fast Fourier transform(FFT) on each row

1D FFT on first row (Note that the absolute value is only displayed and not the actual imaginary number):

Similarly, perform  1D FFT on each row:

NOTE: The figure represents the 1 D FFT on each row and the result is the absolute value of the complex data obtained using FFT.

3.       Perform 1 D Fast Fourier transform on each column.
On the matrix obtained from the previous step, compute 1D FFT column wise.

4.       Display the results obtained.

Flow Chart for Fast Fourier Transform on 2D :

INVERSE FOURIER TRANSFORM:

1.       Perform Inverse Fourier Transform on each column

2.       Perform IFFT on each row

3.       Display the original data

MATLAB CODE:

A=[110 20 140 0 220;

   60 34 23 198 20;

   15 12 126 230 15;

   140 28 10 28 10;

   11 12 19 85 100];

FFT_row = zeros(size(A));

FFT_col = zeros(size(A));

%Perform FFT on each row

for i=1:size(A,1)

FFT_row(i,:) = fft(A(i,:));

end

display(FFT_row);

%display(abs(FFT_row));

%Perform FFT on each column

for i=1:size(A,2)

FFT_col(:,i) = fft(FFT_row(:,i));

end

display(FFT_col);

%display(abs(FFT_col));

%INVERSE FOURIER TRANSFORM

IFFT_row = zeros(size(A));

IFFT_col = zeros(size(A));

%Perform Inverse Fourier Transform on each column

for i=1:size(A,2)

IFFT_col(:,i) = ifft(FFT_col(:,i));

end

%Perform IFFT on each row

for i=1:size(A,2)

IFFT_row(:,i) = ifft(IFFT_col(:,i));

end

display(abs(A))

ALTERNATE METHOD FOR INVERSE FOURIER TRANSFORM:

Instead of using ifft2() or ifft(), we can also use the following method to obtain the original data from the Fast Fourier transformed result :

1.       Obtain the conjugate of the Forward FFT

2.       Perform Forward fast Fourier transform

3.       Obtain the conjugate of the result from step 2.

4.       Divide it by the number of elements present in the matrix

5.       Obtain the original matrix

MATLAB CODE:

Conj_F = conj(F_FFT);

Conj_FFT = fft2(Conj_F);

IFFT_conj = conj(Conj_FFT)/numel(Conj_FFT)

display(abs(IFFT_conj));

Reference: Digital Image Processing  by Rafael C.Gonzalez, fourth Chapter.