Difference between PSD results of the data calculated by the numpy and scipy modules

Question

1. For the calculation of one-dimensional data, I compared the PSD results of numpy and scipy, and found that the calculation results of np.fft.fft() and scipy.signal.periodogram() are the same, but the results of scipy.signal.welch() are not completely different.
2. For one-dimensional time series data, unlike the scipy.signal() module that can directly obtain the PSD value, np.ft.ft() needs to calculate the fast Fourier transform (FFT) of the input sequence and return the complex-valued spectrum, and calculate the spectrum Square the value and divide by a normalization parameter to get the PSD value. However, in this normalization process, I am not sure what parameters to use for normalization. The normalization parameters I see include the total number of samples n, the square of the total number of samples n**2, and the total number of samples multiplied by sampling Frequency n*fs, did not find a too authoritative definition.
3. In addition, the calculation of PSD for two-dimensional matrix data also involves the same normalization parameters of the same dimension data.

Here is the relevant experimental code: scipy

import numpy as np
from scipy.signal import periodogram, welch
import matplotlib.pyplot as plt

# Generate random data
fs = 4000
t = np.arange(0, 3, 1/fs)
y = np.sin(2*np.pi*50*t) + np.cos(2*np.pi*140*t) + np.cos(2*np.pi*300*t) * t ** 2

# Calculate PSD using signal.welch
# f, psd = welch(x, fs=fs, nperseg=256, scaling='density')
f, psd = periodogram(y, fs, scaling='density')

# Plot the PSD versus frequencies
plt.plot(f, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.show()

numpy

# Calculate PSD using np.fft
n = len(y)
y_FFT = np.fft.fft(y)
# psd = np.abs(y_FFT)**2 / (n*fs)
freq = np.fft.fftfreq(n, d=1/fs)
# freq = np.fft.fftshift(freq)
# print(freq)
# Find the indices corresponding to positive frequencies
positive_freq_indices = np.where(freq >= 0)
# Extract the PSD for positive frequencies
psd = np.abs(y_FFT[positive_freq_indices])**2 / (n * fs)
# Plot PSD
plt.plot(freq[positive_freq_indices], psd)
plt.xlabel('Frequency')
plt.ylabel('PSD')
plt.show()

Spectrum analysis based on different PSD results will lead to completely different conclusions. Result drawing visualization The y value changes with time as shown in the figure y value changes

scipy.periodogram scipy.periodogram

scipy.welch scipy.welch

numpy.fft numpy.fft

I'm not sure if the PSD calculation and normalization parameter selection, the calculation of the two-dimensional spatial coordinate axis is correct.

from scipy.fft import fft2, fftfreq

fs = 8000
fft_coefficients = fft2(Two_dim_data)  # Compute the 2D FFT

# Compute the power spectrum by squaring the magnitude of the Fourier coefficients
power_spectrum = np.abs(fft_coefficients) ** 2

# Scale the power spectrum by the appropriate normalization factor
# Normalize by the number of elements in the matrix
two_dim_psd = power_spectrum / (Two_dim_data.shape[0] * fs * Two_dim_data.shape[1] * fs)

# Calculate the spatial frequencies
freq_x = fftfreq(Two_dim_data.shape[0], 1/fs)
freq_y = fftfreq(Two_dim_data.shape[1], 1/fs)

reference_1 reference_2

I wonder why the obtained PSD results vary so much, and what normalization parameter is used when calculating power density from power.