You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: spm_dcm_tfm_image.m
    path of file: ./repos/spm/toolbox/dcm_meeg
    
    the code of the file until where you have to start completion: function spm_dcm_tfm_image(csd,pst,hz,top)
% Display time-frequency complex cross spectra
% FORMAT spm_dcm_tfm_image(csd,top,pst,hz)

    
function spm_dcm_tfm_image(csd,pst,hz,top)
% Display time-frequency complex cross spectra
% FORMAT spm_dcm_tfm_image(csd,top,pst,hz)
% 
% csd - (t x w x n x n): a data array over t time bins, w frequency bins
%                       and n times n channels
% pst - peristimulus time (for plotting)
% Hz  - frequency range (for plotting)
% top - [0/1] switch to display at the top or bottom of the current figure
%__________________________________________________________________________
%
% This routine displays complex cross spectra over peristimulus time as
% images of the absolute values (coherence) and cross covariance functions
% over pairs of channels.
%
% See also: spm_dcm_tfm_responses (and spm_dcm_tfm_transfer) 
%__________________________________________________________________________
 
% Karl Friston
% Copyright (C) 2005-2022 Wellcome Centre for Human Neuroimaging

 
% setup and defaults
%--------------------------------------------------------------------------
two = 2;
if nargin < 2, pst = 1:size(csd,1); end
if nargin < 3, hz  = 1:size(csd,2); end
if nargin < 4, top = 0; two = 1;    end


% plot time frequency responses
%==========================================================================
nc    = size(csd,3);
bands = kron([8; 13; 32],[1 1]);


% evaluate cross covariance function
%--------------------------------------------------------------------------
[ccf,lag] = spm_csd2ccf(csd,hz,1/1024);
lag       = lag*1000;
j         = find(-64 < lag & lag < 64);
lag       = lag(j);
ccf       = ccf(:,j,:,:);


for i = 1:nc    
    
    % spectral power
    %----------------------------------------------------------------------
    for j = i:i
        
        subplot(two*nc,nc,top*nc*nc + (i - 1)*nc + j)
        imagesc(pst,hz,abs(csd(:,:,i,j)).^2');
        title('Spectral density')
        xlabel('peristimulus time (ms)')
        ylabel('Hz'), axis xy
        hold on; plot([pst(1) pst(end