Getting error while running code to get the 'Normalized power spectrum of image'

Am trying to get the normalized power spectrum of an image using matlab, but i get this error whenever i run the program

    ??? Input argument "f" is undefined.
    Error in ==> powerspectrum at 34
    F = fft2(double(f));
The matlab code to the program, please help, am actually a novice when it comes to using matlab, infact i got this code of mathworks website and can't seem to get it to work, any help from anyone would be greatly appreciated, It is part of my final year school project. Someone should please help me out. Thank you

    function [sf nps] = powerspectrum(f, d, ns)
    %%POWERSPECTRUM Get the normalized power spectrum of an image
    %
    %   [SF NPS] = powerspectrum(F) constructs the power spectrum of F by
    %   averaging profiles of the power spectrum image centered at spatial 
    %   frequency 0 mm^-1 and rotated by 1 degree increments. 
    %
    %   [SF NPS] = powerspectrum(F, D) takes as an additional parameter a two
    %   element vector representing pixel size.
    %
    %   [SF NPS] = powerspectrum(F, D, NS) will do the above but take NS spatial
    %   frequency samples.
    %
    %   Example:
    %       f = imread('Lena.tif');
    %       [sf nps] = powerspectrum(f);
    %       semilogy(sf, nps);

    % Pixel dimensions
    if nargin < 2
        d = [1 1];
    end

    % Adjust for nonsquare pixels
    nd = 1./d;
    di = nd/max(nd);

    % Number of samples per line
    if nargin < 3
    ns = 200;
    end

    % Compute the power spectrum
    F = fft2(double(f));
    F = fftshift(F);
    S = abs(F).^2;

    % Coordinates of DC term
    p = size(S)/2;

    % Sampling matrix
    PS = zeros(360, ns);

    % Sample along lines rotated about DC term
    for theta = 1:360
    u = [cosd(theta) sind(theta)];
    q = p.*(1 + u).*di;
    PS(theta,:) = linesample(S, p, q, ns);
    end

    % Average and normalize
    ps = mean(PS,1);
    nps = ps/max(ps);

    % Spatial frequencies
    sfu = 1/(2*min(d)); % Nyquist criteria
    sfi = sfu/ns;
    sf = 0:sfi:sfu-sfi;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function v = linesample(f, p, q, n, method)
    %%LINESAMPLE Sample image intensities along a line 
    %
    %   V = linesample(F, P, Q) will sample the image F between points P and Q. 
    %
    %   V = linesample(F, P, Q, N) will do the above but take N number of evenly
    %   spaced samples. 
    %
    %   V = linesample(..., METHOD) will do the above but use the specified
    %   interpolation method. Bilinear (the default) and nearest neighbor
    %   interpolation are currently supported.

    % Calc number of samples, if necessary
    if nargin < 4
    n = round(sqrt(sum((p-q).^2)));
    end

    % Calculate step size
    step = (q-p)/n;

    % Get points at which to sample
    for d = 1:2 
    if p(d) == q(d)
        x{d} = p(d)*ones(1,n);
    else
        x{d} = p(d):step(d):q(d);
      end
    end

    % Sample 
    if nargin < 5 || strcmp(method, 'bilinear')
      v = bilinear(f, x{1}(1:n), x{2}(1:n));
    else 
      v = nearestneighbor(f, x{1}(1:n), x{2}(1:n));
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function zi = bilinear(f, xi, yi)
    %%Bilinear interpolation

    d = size(f);

    xu = ceil(xi);
    xl = floor(xi);

    yu = ceil(yi);
    yl = floor(yi);

    % Fix out of range subscripts
    xu(xu > d(1)) = d(1);
    xl(xl == 0) = 1;
    yu(yu > d(2)) = d(2);
    yl(yl == 0) = 1;

    xuyu = sub2ind(d, xu, yu);
    xlyu = sub2ind(d, xl, yu);
    xuyl = sub2ind(d, xu, yl);
    xlyl = sub2ind(d, xl, yl);

    a = (xi-xl).*(f(xuyu)-f(xlyu))+f(xlyu);
    b = (xi-xl).*(f(xuyl)-f(xlyl))+f(xlyl);

    zi = (yi-yl).*(a-b)+b;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function zi = nearestneighbor(f, xi, yi)
    %%Nearest neighbor interpolation

    d = size(f);

    x = round(xi);
    y = round(yi);

    x(x > d(1)) = d(1);
    x(x == 0) = 1;
    y(y > d(2)) = d(2);
    y(y == 0) = 1;

    idx = sub2ind(size(f), x, y);
    zi = f(idx);