I am trying to take fft and dtft of a grey scale image from first principals(so without using built it fft or dtft functions) using python

I need to write a function that computes the fft and another to compute the dft of a a composite sinusoidal signal and to then plot the phase and magnitude spectrum there of. I then need to use the same functions to compute the fft and dft of a greyscale image. I have the functions and can get the transforms of a composite sine wave but I am not sure how to apply this to an image. I believe that the pixels in an image are stored as a 2D array where as the sine function is only 1D

here is my code:

import numpy as ny 
import matplotlib.pyplot as mplot 
import time as tm # used to claculate time taken to perform functions 

# Variables
sample_rate = 0
sample_period = 0
time = 0 # Total time
T = 0
n = 0
freq = 0
y = 0 # Used to build signal to apply the DTF to
Y = 0 # Signal after the DTF
#***

# Main
sample_rate = 128
sample_period = 1/sample_rate
time = ny.arange(0,2,sample_period) # Sample set up to got from 0 - 2 seconds

# Signal 1
freq = 1
y = ny.sin(2*ny.pi*freq*time)
# Signal 2
freq = 3
y += 2*ny.sin(2*ny.pi*freq*time)
# Signal 3
freq = 5
y += 3*ny.sin(2*ny.pi*freq*time)
# Signal 4
freq = 7
y += 0.5*ny.sin(2*ny.pi*freq*time)

#time domain plot of signal
mplot.figure(figsize = (6,4))
mplot.plot(time, y, 'b')
mplot.title('Plot of signal y')
mplot.xlabel('Time (s)')
mplot.ylabel('Amplitude')
#mplot.show()
#***

# Functions
# Discrete Fourier Transform
def DFT(y):
    n = ny.arange(len(y))
    k = n.reshape((len(y), 1))
    e = ny.exp(-2j * ny.pi * k * n / len(y))
    
    Y = ny.dot(e,y)
    
    return Y

# Fast Fourier Transform
def FFT(y):

    length = len(y)
    if length == 1:
        return y
    else:
        even = FFT(y[::2])
        odd = FFT(y[1::2])
        fact = ny.exp(-2j * ny.pi * ny.arange(length) / length)
        
        Y = ny.concatenate(( even + odd*fact[:int(length/2)] , even + odd*fact[int(length/2):] )) # Joining the arrays together
        return Y

#***
#get DFT of composite signal
dftstart = tm.time()
Y = DFT(y)
dftend = tm.time()
dfttime = dftend - dftstart
# calculate the frequency
n = ny.arange(len(Y))
T = len(Y)/sample_rate
freq = n/T 

#plot magnitude spctrum of DTFT of composite signal 
mplot.figure(figsize = (6,4))
mplot.stem(freq, abs(Y), 'r')
mplot.title('Discrete Fourier Transform of signal y')
mplot.xlabel('Frequency (Hz)')
mplot.ylabel('DFT Amplitude |Y(e^jw)|')
#mplot.show()

#plot phase response of DFT
mplot.figure(figsize = (6,4))
mplot.plot(freq, Y/abs(Y), 'r')
mplot.title('phase spectrum of DFT of signal y')
mplot.xlabel('Frequency (Hz)')
mplot.ylabel('Phase(rad/s)')


#get fft of composite signal
fftstart = tm.time()
Y = FFT(y)
fftend = tm.time()
ffttime = fftend - fftstart
# calculate the frequency
length = len(Y)
n = ny.arange(length)
T = len(Y)/sample_rate
freq = n/T 

#plot magnitude spectrum of fft of composite signal 
mplot.figure(figsize = (6,4))
mplot.stem(freq, abs(Y), 'r')
mplot.title('Fast Fourier Transform of signal y')
mplot.xlabel('Frequency (Hz)')
mplot.ylabel('FFT Amplitude |X(e^jw)|')
#mplot.show()

#plot phase response of fft
mplot.figure(figsize = (6,4))
mplot.plot(freq, Y/abs(Y), 'b')
mplot.title('phase response of FFT of signal y')
mplot.xlabel('Frequency (Hz)')
mplot.ylabel('Phase(rad/s)')
mplot.show()

print("time to execute dft: " + str(dfttime))
print("time to execute fft: " + str(ffttime))

If someone could advise me on how to apply this to a greyscale image it would be greatly appreciated