Download Homework 3 Solutions - Digital Image Processing | ECE 533 and more Assignments Digital Signal Processing in PDF only on Docsity!
Digital Image Processing ECE 533
Solutions Assignment 3
Date: March 11, 2008
Department of Electrical and Computing Engineering, University of New Mexico.
Professor Majeed Hayat, [email protected]
March 11, 2008
1 The sampling theorem
Part 1:
- The minimum sampling rate necessary to assure perfect reconstruction, i.e., the Nyquist rate, is equal to twice the highest frequency found in the signal spectrum. In the case of y(t) = cos(2Ļt), all the spectrum is concentrated at the frequencies ā 2 Ļ and 2Ļ, thus the minimum sampling rate is 4Ļ rad/s, or 2 Hertz.
- For a sampling rate fs = T ā^1 greater than the Nyquist rate, the spectrum (Fourier transform) of the sampled signal ys(t) = āā n=āā y(nT )Ī“(t ā nT ) is composed of a series of copies of the spectrum of y(t) multiplied by fs and centered in 2Ļnfs, n =... , ā 1 , 0 , 1 ,.. .. Considering the signal y(t) = cos(2Ļt), the reconstruction filter has to filter in the copy that is around 0, while filtering out the other copies. For the signal in question, the copy of the spectrum centered in 2Ļfs has a component in 2Ļ(fs ā 1), thus B has to be in the range 2Ļ ā¤ B < 2 Ļ(fs ā 1). If B is less
than this lower limit, then the reconstructed signal would be zero, if it is beyond the upper limit, other frequency components that didnāt exist in the original signal will appear in the reconstructed signal.
- From the discussion in the answer to the previous question, perfect reconstruction is assured if 1 ā¤ p < T ā^1 ā 1.
- We have the following issues
Aliasing: it doesnāt occur in none of the cases, since the sampling rate is fs = T ā^1 = 52 > 2. Suitability of the bandwidth of the filter: from the previous discussion, 1 ā¤ p < T ā^1 ā 1 = 32. Hence, the bandwidth is suitable only in the cases p = 65 , 1. Discrepancies in the amplitude: the reconstructed signal should show discrep- ancies when p is outside of the range 1 ā¤ p < 32. Effects near the ends of the time interval [ā 2 , 2]: Except for the sampling in- stants, at any instant of time the reconstructed signal receives contributions from all the samples. This happens because the sinc functions that compose the reconstructed signal are non-zero almost everywhere. Hence, to have a per- fectly reconstructed signal it would be necessary to take into account all the infinite samples. However, the further away a sample is from the time-instant where the signal is being reconstructed, the lesser contribution it represents. That is why the points closer to the ends of the time interval may present greater discrepancy between the original and reconstructed signal. Frequency components of the reconstructed signal: When p < 1, all the fre- quency components of the sampled signal will be filtered out, so the resulting reconstructed signal should have no components at all. When 1 ā¤ p < 32 , only the desired frequency component should be present. When p = 32 , there should be 2 frequency components: at 1.0 and 32 Hertz. When p = 2 (respectively, p = 4), there should be frequency components at 1.0 and 32 (respectively, 1.0, (^32) and 92 ) Hz.
ā4^ ā
0 2
4 ā2 ā
(^20) 40
v u
|F(u,v)|
Figure 1: Left: Sampled and reconstructed versions of z(x, y) = cos(2Ļ(x + y)) for fs,x = fs,x = 52 cmā^1 and px = py = 65. Right: The magnitude of the 2D Fourier transform of z(x, y). .
2 Fourier transform, orthogonal representation of im-
ages and simple image compression
In order to derive a theoretical expression to calculate the coefficients Ī»m,n, we start with the equation I = UĪVT^ and perform the following matrix operation
I = UĪVT^ (5) UT^ IV = UT^ UĪVT^ V (6) UT^ IV = Ī, (7)
where the last equation holds because we know that for unitary transformations: UUT^ = I and VVT^ = I. For the MATLAB implementation see the script provided below. In the first compression algorithm we simply discard N high-frequency components in each spatial-frequency direction. In Fig. 2, some sample images are given. For the first image, football.jpg, N was set to 90 and for the second image, cameraman.tif, N = 80. These values for N mean that all those frequencies above the N th FFT point are defined to be zero. For this kind of compression algorithm we achieve compression because we
reduce the number of frequency components we have to store or transmit to the decoding side. Note the smoothing effect of the compression algorithm on the images. The second algorithm keeps all the spatial-frequency components but truncates the FFT coefficients using no more than M decimal places per coefficient. For the images shown in Fig. 2, M = 1. In this second case, we achieve compression because instead of using floating point numbers for representing the FFT coefficients, we can multiply such numbers by 10 M^ and use an integer representation. Note that there is no smoothing effect on the images (why?), but some granular noise appears on top of the images.
% ECE-533, Fall 2008 % MATLAB program for computing the orthogonal representation of an image clear all; clc Image2Use=2; switch Image2Use case 1 Image=imread(āpout.tifā); Filename=āpoutā; % Number of frequencies to eliminate Freq2Eliminate=90; case 2 Image=imread(āfootball.jpgā); Image=rgb2gray(Image); Filename=āfootballā; % Number of frequencies to eliminate Freq2Eliminate=80; case 3 Image=imread(ātrees.tifā); Filename=ātreesā; % Number of frequencies to eliminate Freq2Eliminate=50; case 4
Image=imread(āsaturn.pngā); Image=rgb2gray(Image); Filename=āsaturnā; % Number of frequencies to eliminate Freq2Eliminate=80; case 5 Image=imread(ācameraman.tifā); Filename=ācameramanā; % Number of frequencies to eliminate Freq2Eliminate=80; end % Recast the image as a double object Image=double(Image); % Compute the sizes M=size(Image,1); N=size(Image,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Computation of the linear components of the image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Define the matrices U and V [x,y]=meshgrid(0:(M-1),0:(M-1)); U=mod(x.y,M); U=exp(2pii/MU)/sqrt(M); [x,y]=meshgrid(0:(N-1),0:(N-1)); V=mod(x.y,N); V=exp(2pii/NV)/sqrt(N); % Compute the coefficients lambda Lambda=UāImageV; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Compute the 2D FFT of the image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the 2D FFT of the image IMAGE=fft2(Image)/(MN); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPRESSION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Compress the image dropping some terms IMAGEComp=IMAGE; IMAGEComp=fftshift(IMAGEComp); % Eliminate some high frequency components IMAGEComp(1:Freq2Eliminate,:)=0; IMAGEComp((end-Freq2Eliminate):end,:)=0; IMAGEComp(:,1:Freq2Eliminate)=0; IMAGEComp(:,(end-Freq2Eliminate-1):end)=0; ImageComp=ifft2(fftshift(IMAGEComp)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Compress the image truncating all some terms Decs=1; % Number of decimal positions to keep IMAGEComp2=round(IMAGE10^Decs)/10^Decs; ImageComp2=ifft2(IMAGEComp2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Show images figure(1); subplot(231) imshow(Image,[]); xlabel(āOriginal Imageā,āFontSizeā,14); subplot(232) imshow(log10(abs(Lambda+eps.^2)),[]); xlabel(āFourier decompositionā,āFontSizeā,14); subplot(233) imshow(log10(abs(IMAGE+eps).^2),[]);
