Hi ZW,
I did slight modification to your code, it is working fine for me. Converting "psd_range" to decibels (dB) for better visualization.
Let me know if there is anything specific you want to know?
% Parameter settings
Fs = 1e6; % Sampling rate 1 MHz
t_end = 1; % Signal duration 1 second
t = 0:1/Fs:t_end-1/Fs; % Time vector
% Frequencies
frequencies = [1e3, 2e3, 5e3, 13e3]; % Frequencies of the sine waves
amplitude = 1; % Amplitude of the sine waves
% Generate composite signal
combined_signal = sum(amplitude * sin(2 * pi * frequencies' * t), 1);
% PSD power spectrum calculation
[psd, f] = pwelch(combined_signal, [], [], 4096*4096, Fs); % Calculate power spectrum density
% Limit the frequency range to 20-20000 Hz
f_range = f >= 20 & f <= 20000;
psd_range = psd(f_range);
f_range = f(f_range);
% Plotting
figure;
plot(f_range, 10*log10(psd_range)); % Convert to dB
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
title('Power Spectral Density');
grid on;
% Calculate signal power
signal_power = sum(abs(combined_signal).^2) / length(combined_signal);
P = Fs / 2 / length(psd_range) * (sum(psd_range) - 0.5 * (psd_range(1) + psd_range(end)));
