%%: Edited Diary of Matlab Demo Session, in ECEU692, Intro to SSI Class %%: Demo presented in class #7, 2 Feb, 2004 %%: %%: All comments added class, and all are marked by a preceeding '%' %%: %%: D.H. Brooks %%: v1.0 %%: Feb 2004 %% 'what' command shows all m-files and .mat files in current directory: what M-files in the current directory C:\document\working\10471\dana lowpass_image space_cos %%: First we look at spatial cosines to illustrate frequency ideas in space %%: and in the Fourier domain. more on % makes Matlab pause between screen-fuls of text for a 'return' %% display the m-file (also given in Lecture 6 Slide which creates four 2D %% signals which are separable products of 1D discrete cosines %% all 4 use a vector outer product (column vector by a row vector) to %% create the matrix that represents the 2D signal %% Note the use of the 'ones()' function to get a matrix or vector of all %% ones. Also note the use of the (:) after a variable to force it to be a %% column vector, and the .' to get the transponse. type space_cos % 'type' types the content of an m-file to the screen x = cos(2*pi/32*[0:127]); y = cos(2*pi/16*[0:127]); z = cos(2*pi/4*[0:127]); horcos = ones(128,1)*(x(:).'); vercos = y(:)*ones(1,128); horver = y(:) * (x(:).'); horfver = z(:) * (x(:).'); % run space_cos.m space_cos % who shows all variables in workspace. whos shows a full description of % all variables. who Your variables are: horcos horver x z horfver vercos y whos Name Size Bytes Class horcos 128x128 131072 double array horfver 128x128 131072 double array horver 128x128 131072 double array vercos 128x128 131072 double array x 1x128 1024 double array y 1x128 1024 double array z 1x128 1024 double array Grand total is 65920 elements using 527360 bytes % display the results using 'mesh' mesh(horcos) mesh(horfver) % to examine the frequency content, use fft2 to get a 2D Fourier Transform % (technically, the DFT, computes using an fft algorithm) fhorcos = fft2(horcos); whos Name Size Bytes Class ans 1x1 8 double array fhorcos 128x128 262144 double array (complex) horcos 128x128 131072 double array horfver 128x128 131072 double array horver 128x128 131072 double array vercos 128x128 131072 double array x 1x128 1024 double array y 1x128 1024 double array z 1x128 1024 double array Grand total is 82305 elements using 789512 bytes % Note that since the DFT is complex we need to use 'abs()' (or 'angle()', % 'real()', or 'imag()') for display purposes mesh(abs(fhorcos)) % Note that the DFT runs from 0 to pi (in both dimensions for 2D signals) % in frequency. To see the DC (0 frequency) value in the middle we use % fftshift. mesh(abs(fftshift(fhorcos))) fvercos = fft2(vercos); mesh(abs(fftshift(fvercos))) fhorver = fft2(horver); mesh(abs(fftshift(fhorver))) fhorver = fft2(horfver); fhorfver = fft2(horfver); mesh(abs(fftshift(fhorfver))) %%: Now we look at the frequency content of a 'real-world' image (our old %%: friend the California band) and what happens when we modify it % ls shows full contents of current directory ls . calband.tif lowpass_image.m .. class_020204.txt space_cos.m % read in image into variable 'calband' [calband,b] = imread('calband.tif', 'tif'); whos Name Size Bytes Class ans 1x1 8 double array b 256x3 6144 double array calband 696x608 423168 uint8 array fhorcos 128x128 262144 double array (complex) fhorfver 128x128 262144 double array (complex) fhorver 128x128 262144 double array (complex) fvercos 128x128 262144 double array (complex) horcos 128x128 131072 double array horfver 128x128 131072 double array horver 128x128 131072 double array vercos 128x128 131072 double array x 1x128 1024 double array y 1x128 1024 double array z 1x128 1024 double array Grand total is 555393 elements using 2005256 bytes imagesc(calband) axis('image'),colormap('gray') % note that using axis('image') causes the x and y axis to have the same % length scale in the display % let's look at a function I wrote to put zeros at a specified number of % frequencies: type lowpass_image function smoother_transform = lowpass_image(input_transform,cutoff); [xdim,ydim] =size(input_transform); smoother_transform = input_transform; smoother_transform(1:cutoff(1),:) = 0; smoother_transform(:,1:cutoff(2)) = 0; smoother_transform(xdim - cutoff(1)+1:xdim,:) = 0; smoother_transform(:,ydim - cutoff(2)+1:ydim) = 0; % Note in particular the use of the word 'function' as the first word in % the m-file to tell Matlab that this is a function (with arguments passed % in and out) and not a script. Note also the use of the : operator here to % do 2 things: denote a range of integers (as in 1:cutoff(1)) and denote % the entire range of rows or columns in a matrix (as in the second colon % in smoother_transform(1:cutoff(1),:) = 0; which means take all columns of % the specified rows. % Since this function is designed to operate as a crude low-pass filter % when applied to the Fourier Transform of an image after it is arranged % with the DC value at the center, first we take the 2D-DFT of the calband % image fcalband = fft2(calband); who Your variables are: ans fhorcos horcos x b fhorfver horfver y calband fhorver horver z fcalband fvercos vercos whos Name Size Bytes Class ans 1x1 8 double array b 256x3 6144 double array calband 696x608 423168 uint8 array fcalband 696x608 6770688 double array (complex) fhorcos 128x128 262144 double array (complex) fhorfver 128x128 262144 double array (complex) fhorver 128x128 262144 double array (complex) fvercos 128x128 262144 double array (complex) horcos 128x128 131072 double array horfver 128x128 131072 double array horver 128x128 131072 double array vercos 128x128 131072 double array x 1x128 1024 double array y 1x128 1024 double array z 1x128 1024 double array Grand total is 978561 elements using 8775944 bytes % Now we explore some ways to display this Fourier Transform: imagesc(abs(fcalband)) mesh(abs(fcalband)) mesh(fftshift(abs(fcalband))) mesh(log10(fftshift(abs(fcalband)))) % now we display the effects of the lowpass_image function on a small example lowpass_image(ones(10,8),[2,1]) ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 % and now we apply it to the fftshift'ed 2D DFT of calband fsmoothed_band = lowpass_image(fftshift(fcalband),[100 100]); smoothed_band = ifft2(ifftshift(fsmoothed_band)); whos Name Size Bytes Class ans 10x8 640 double array b 256x3 6144 double array calband 696x608 423168 uint8 array fcalband 696x608 6770688 double array (complex) fhorcos 128x128 262144 double array (complex) fhorfver 128x128 262144 double array (complex) fhorver 128x128 262144 double array (complex) fsmoothed_band 696x608 6770688 double array (complex) fvercos 128x128 262144 double array (complex) horcos 128x128 131072 double array horfver 128x128 131072 double array horver 128x128 131072 double array vercos 128x128 131072 double array x 1x128 1024 double array y 1x128 1024 double array z 1x128 1024 double array Grand total is 1401808 elements using 15547264 bytes --more-- whos Name Size Bytes Class ans 10x8 640 double array b 256x3 6144 double array calband 696x608 423168 uint8 array fcalband 696x608 6770688 double array (complex) fhorcos 128x128 262144 double array (complex) fhorfver 128x128 262144 double array (complex) fhorver 128x128 262144 double array (complex) fsmoothed_band 696x608 6770688 double array (complex) fvercos 128x128 262144 double array (complex) horcos 128x128 131072 double array horfver 128x128 131072 double array horver 128x128 131072 double array smoothed_band 696x608 6770688 double array (complex) vercos 128x128 131072 double array x 1x128 1024 double array y 1x128 1024 double array z 1x128 1024 double array Grand total is 1824976 elements using 22317952 bytes % since we have mucked around with the 2D DFT rather carelessly, not that % the inverse DFT gives a complex-valued rather than a real-valued % result. One way to see how strong this effect is, is to compare the max % of the imaginary part to the max of the real part. Note the use of the % nested 'max()' functions since max() of an matrix in Matlab returns a % vector containing the max of each column. max(max(imag(smoothed_band))) ans = 3.0221 max(max(real(smoothed_band))) ans = 309.1659 % Satisfied that the real part is 2 orders of magnitdue (at least its % maximum is) greater than the imaginary part, we proceed using the abs() % function to display real-valued quantities imagesc(abs(smoothed_band)) axis('image'),colormap('gray') figure,imagesc(calband) axis('image'),colormap('gray') % now let's try a more dramatic amount of filtering fsmoothed_band_31 = lowpass_image(fftshift(fcalband),[300 100]); smoothed_band_31 = ifft2(ifftshift(fsmoothed_band_31)); imagesc(abs(smoothed_band_31)) axis('image'),colormap('gray') diary off